MPS & DMRG Hands-On Session ASC Summer School 2017 at LMU Munich

Table of Contents

1 Introductory remarks

These notes are meant to go alongside the hands-on session and to facilitate looking up a command which already scrolled off the screen, to re-read at a later time and/or to follow along individually at a faster or slower pace. The two sections marked [extra] likely won’t be discussed in the main section due to limited time. The first two parts will be discussed on Monday, the fourth and fifth part on Tuesday.

Commands to be entered will be typeset like this:

1: $ command-to-be-entered argument1 argument2

and their output will be given like this:

2: some output
3: some more output

Line numbers are continuous through the document with the exception of boxes meant to be copy-pasted.

2 Basics of matrix-product operators & states

2.1 Files and executables

The SyTen toolkit primarily knows two kinds of files: lattice files, which represent a concrete, finite physical system including its operators, local bases etc. and state files, which store the tensor-product representation of a single, pure quantum state.

As such, the first step in any calculation is always to generate a lattice file. This can be done using one of the lattice generators, executables which, when run, generate a lattice file for a physical system of a certain size, interactions, local bases etc. Afterwards, generic operators may read in this lattice file to e.g. construct a random state compatible with the physical system, print tensor product operator representations stored in the lattice, calculate expectation values etc. We will go over these step by step.

2.2 Generating a lattice representation & matrix-product operators

While it is relatively easy to write a lattice generator using the existing library, we will only use pre-existing ones here. One of those is the syten-sql-mps-spin executable. Like all executables, it will offer a detailed help when passed the --help flag:

4: $ syten-sql-mps-spin --help
 5: Spin square lattice.
 6: Defaults to creating a 10×1 square lattice in the file spin-sql-<L>x<W>-<S>-<sym>.lat.
 7: Usage: syten-sql-mps-spin [options]
 8: Allowed options:
 9:   -l [ --length ] arg (=10)         length of the system (long direction)
10:   -w [ --width ] arg (=1)           width of the system (short direction)
11:   -z [ --sz ] arg (=0.5)            spin size on each site
12: […snip…]

Using this tool, we can generate a lattice of four \(\mathrm{U}(1)\)-symmetric \(S=\frac{1}{2}\) spins in file sp4.lat:

13: $ syten-sql-mps-spin --sym u1 -l 4 -z 0.5 -o sp4.lat
14: I: Creating a 4×1 spin lattice in file sp4.lat, using mapping variant Z and shift 0 with u1 symmetries.
15:   Ordering is defined by new[old]: [0, 1, 2, 3].
16: ....
17: 2017.08.09T15:42:29.192 3.51510e-02 I Working on operator Hp
18: 2017.08.09T15:42:29.192 3.54550e-02 I Working on operator Hp_xp
19: 2017.08.09T15:42:29.192 3.54870e-02 I Working on operator Hz
20: 2017.08.09T15:42:29.192 3.56560e-02 I Working on operator Hz_s
21: 2017.08.09T15:42:29.199 4.01110e-02 I Working on operator Hz_xp
22: Saving lattice…
23: Saved lattice into file sp4.lat.

This file now represents our lattice of four spin one-halfs and we can inspect the stored information using syten-info as well as syten-print. First, some general information on the lattice:

24: $ syten-info sp4.lat
25: 2017.08.09T15:44:52.013 2.27260e-02 I Loaded MP lattice
26: == MPS Lattice Information:
27:    U(1)-symmetric S=5.000e-01 spin chain on a square lattice of dimension 4×1 with mapping:
28:           (3)
29:    (  0)   0  (  0)
30:    (  1)   1  (  1)
31:    (  2)   2  (  2)
32:    (  3)   3  (  3)
33:            (0)
34: 
35:    Vacuum basis:         [<1, [U(1)_z:{+0}@1]>]
36:    Number of sites:      4
37:    Uniform lattice:      Yes
38:    Local site basis:     [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
39: -- Local site operators:
40:    sm                    s^-_i
41:    sp                    s^+_i
42:    sz                    s^z_i
43: -- Global operators:
44:    Hp                    \sum_{<ij>} 1/2 [S^+_i S^-_j + S^-_i S^+_j]   uniform coupling in x- & y-direction [obc]
45:    Hp_xp                 \sum_{<ij>} 1/2 [S^+_i S^-_j + S^-_i S^+_j]   uniform coupling in x- & y-direction [add. term for length-wise pbc]
46:    Hz                    \sum_{<ij>} S^z_i S^z_j                       uniform coupling in z-direction [obc]
47:    Hz_s                  \sum_{i} (-1)^i S^z_i                         staggered magnetic field terms in z-direction
48:    Hz_xp                 \sum_{<ij>} S^z_i S^z_j                       uniform coupling in z-direction [add. term for length-wise pbc]
49:    I                     global identity operator
50:    Proj_Down             \prod_{i=1}^N (0.5 - S^z_i)
51:    Proj_Up               \prod_{i=1}^N (0.5 + S^z_i)
52: == End.

And on specific operators:

53: $ syten-info sp4.lat:Hp
54: 2017.08.09T15:45:39.858 2.39840e-02 I Loaded MP operator.
55: 2017.08.09T15:45:39.859 2.41400e-02 I MP-Operator on lattice of 4 sites.
56: 2017.08.09T15:45:39.859 2.41960e-02 I Vacuum basis (right edge): [<1, [U(1)_z:{+0}@1]>]
57: 2017.08.09T15:45:39.859 2.42600e-02 I Sector basis (left edge):  [<1, [U(1)_z:{+0}@1]>]
58: 2017.08.09T15:45:39.859 2.42820e-02 I Physical basis on site 1:  [<1, [U(1)_z:{+0.5}@1]>, <1, [U(1)_z:{-0.5}@1]>]
59: $ syten-info sp4.lat:'sz 0 @'
60: 2017.08.09T15:45:43.745 2.55920e-02 I Loaded MP operator.
61: 2017.08.09T15:45:43.745 2.57560e-02 I MP-Operator on lattice of 4 sites.
62: 2017.08.09T15:45:43.746 2.58040e-02 I Vacuum basis (right edge): [<1, [U(1)_z:{+0}@1]>]
63: 2017.08.09T15:45:43.746 2.58660e-02 I Sector basis (left edge):  [<1, [U(1)_z:{+0}@1]>]
64: 2017.08.09T15:45:43.746 2.59050e-02 I Physical basis on site 1:  [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
65: $ syten-info sp4.lat:'sp 0 @'
66: 2017.08.09T15:45:45.799 2.32210e-02 I Loaded MP operator.
67: 2017.08.09T15:45:45.799 2.34410e-02 I MP-Operator on lattice of 4 sites.
68: 2017.08.09T15:45:45.799 2.34830e-02 I Vacuum basis (right edge): [<1, [U(1)_z:{+0}@1]>]
69: 2017.08.09T15:45:45.799 2.35290e-02 I Sector basis (left edge):  [<1, [U(1)_z:{+1}@1]>]
70: 2017.08.09T15:45:45.799 2.35520e-02 I Physical basis on site 1:  [<1, [U(1)_z:{+0.5}@1]>]

We can also print information about the dimensions of each MPO block using -b. Columns are the bond, number of quantum number sectors, effective dimension, total dimension, average dense block size and the size of the largest dense block.

71: $ syten-info -qb sp4.lat:'Hp'
72: (0,1)         3         4         4         1         2
73: (1,2)         3         4         4         1         2
74: (2,3)         3         4         4         1         2
75: $ syten-info -qb sp4.lat:'Hz'
76: (0,1)         1         3         3         3         3
77: (1,2)         1         3         3         3         3
78: (2,3)         1         3         3         3         3
79: $ syten-info -qb sp4.lat:'Hz Hp +'
80: (0,1)         3         4         4         1         2
81: (1,2)         3         5         5         2         3
82: (2,3)         3         4         4         1         2

Printing the representation of a specific operator is also possible, though potentially tedious (with comments further down and closer to the location they reference):

83: $ syten-print sp4.lat:'sz 1 @'
84: Loaded MP operator sp4.lat:sz 1 @.
85: Site: 0
86: ========================================
87: Tensor of rank 4 with the following bases:
88: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
89: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
90: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
91: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
92: Composed of the following 2 blocks:

The ‘header’ of the tensor on the first site. Leg 1 is the RHS MPO leg, leg 2 the LHS MPO leg. Leg 3 is the upper physical leg and leg 4 the lower physical leg. A MPS would connect to the fourth leg. Since this MPO transforms trivially, the two MPO legs only carry the \(S^z = 0\) quantum number sector.

93: ----------------------------------------
94: Tensor block of rank 4 and 1 symmetries. isZero = +0.
95: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
96: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1

The first block maps the \(S^z=-\frac{1}{2}\) to \(S^z=-\frac{1}{2}\) with prefactor 1.

 97: ----------------------------------------
 98: Tensor block of rank 4 and 1 symmetries. isZero = +0.
 99: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
100: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1

The second block does the same in the \(S^z=\frac{1}{2}\) sector, this is an identity MPO component.

101: ========================================
102: 
103: Site: 1
104: ========================================
105: Tensor of rank 4 with the following bases:
106: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
107: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
108: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
109: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
110: Composed of the following 2 blocks:
111: ----------------------------------------
112: Tensor block of rank 4 and 1 symmetries. isZero = +0.
113: Reduced tensor is: (-5.0000e-01,+0.0000e+00)
114: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
115: ----------------------------------------
116: Tensor block of rank 4 and 1 symmetries. isZero = +0.
117: Reduced tensor is: (+5.0000e-01,+0.0000e+00)
118: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1

This site is more interesting, as the \(\hat S^z_1\) actually acts here – the blocks gain a prefactor \(\pm \frac{1}{2}\) depending on their quantum number sectors.

119: ========================================
120: 
121: Site: 2
122: ========================================
123: Tensor of rank 4 with the following bases:
124: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
125: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
126: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
127: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
128: Composed of the following 2 blocks:
129: ----------------------------------------
130: Tensor block of rank 4 and 1 symmetries. isZero = +0.
131: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
132: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
133: ----------------------------------------
134: Tensor block of rank 4 and 1 symmetries. isZero = +0.
135: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
136: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
137: ========================================
138: 
139: Site: 3
140: ========================================
141: Tensor of rank 4 with the following bases:
142: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
143: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
144: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
145: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
146: Composed of the following 2 blocks:
147: ----------------------------------------
148: Tensor block of rank 4 and 1 symmetries. isZero = +0.
149: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
150: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
151: ----------------------------------------
152: Tensor block of rank 4 and 1 symmetries. isZero = +0.
153: Reduced tensor is: (+1.0000e+00,+0.0000e+00)
154: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
155: ========================================

The remaining two sectors are just identities once more.

syten-print does not compress its arguments by default, so we have to apply the additional _truncate operation if our MPO is more complicated:

156: $ syten-print sp4.lat:'Hp Hz + _truncate'
157: Loaded MP operator sp4.lat:Hp Hz + _truncate.
158: Site: 0
159: ========================================
160: Tensor of rank 4 with the following bases:
161: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <2, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
162: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
163: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
164: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
165: Composed of the following 4 blocks:
166: ----------------------------------------
167: Tensor block of rank 4 and 1 symmetries. isZero = +0.
168: Reduced tensor is:
169: Start: {2, 1, 1, 1}/2
170: [0, 0, r, c]
171: 	(-5.0000e-01,+0.0000e+00)	
172: [1, 0, r, c]
173: 	(+1.0000e+00,+0.0000e+00)	
174: End: {2, 1, 1, 1}/2
175: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
176: ----------------------------------------
177: Tensor block of rank 4 and 1 symmetries. isZero = +0.
178: Reduced tensor is:
179: Start: {2, 1, 1, 1}/2
180: [0, 0, r, c]
181: 	(+5.0000e-01,+0.0000e+00)	
182: [1, 0, r, c]
183: 	(+1.0000e+00,+0.0000e+00)	
184: End: {2, 1, 1, 1}/2
185: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
186: ----------------------------------------
187: Tensor block of rank 4 and 1 symmetries. isZero = +0.
188: Reduced tensor is: (+7.0711e-01,+0.0000e+00)
189: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{-1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
190: ----------------------------------------
191: Tensor block of rank 4 and 1 symmetries. isZero = +0.
192: Reduced tensor is: (+7.0711e-01,+0.0000e+00)
193: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
194: ========================================

This looks different from the usual representation, because the tensor legs are intertwined in a funny way. However, when disentangling it it’s precisely [ S^z ; 1 ; (1/sqrt(2) S^+) ; (1/sqrt(2) S^-) ]

195: Site: 1
196: ========================================
197: Tensor of rank 4 with the following bases:
198: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <3, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
199: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <2, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
200: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
201: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
202: Composed of the following 6 blocks:
203: ----------------------------------------
204: Tensor block of rank 4 and 1 symmetries. isZero = +0.
205: Reduced tensor is:
206: Start: {3, 2, 1, 1}/6
207: [0, 0, r, c]
208: 	(-5.0000e-01,+0.0000e+00) A	
209: [0, 1, r, c]
210: 	(+0.0000e+00,+0.0000e+00) B	
211: [1, 0, r, c]
212: 	(+0.0000e+00,+0.0000e+00) C	
213: [1, 1, r, c]
214: 	(-5.0000e-01,+0.0000e+00) D	
215: [2, 0, r, c]
216: 	(+0.0000e+00,+0.0000e+00) E	
217: [2, 1, r, c]
218: 	(+1.0000e+00,+0.0000e+00) F	
219: End: {3, 2, 1, 1}/6
220: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
221: ----------------------------------------
222: Tensor block of rank 4 and 1 symmetries. isZero = +0.
223: Reduced tensor is:
224: Start: {3, 2, 1, 1}/6
225: [0, 0, r, c]
226: 	(+5.0000e-01,+0.0000e+00) G	
227: [0, 1, r, c]
228: 	(+0.0000e+00,+0.0000e+00) H	
229: [1, 0, r, c]
230: 	(+0.0000e+00,+0.0000e+00) I	
231: [1, 1, r, c]
232: 	(+5.0000e-01,+0.0000e+00) J	
233: [2, 0, r, c]
234: 	(+0.0000e+00,+0.0000e+00) K	
235: [2, 1, r, c]
236: 	(+1.0000e+00,+0.0000e+00) L	
237: End: {3, 2, 1, 1}/6
238: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
239: ----------------------------------------
240: Tensor block of rank 4 and 1 symmetries. isZero = +0.
241: Reduced tensor is:
242: Start: {1, 2, 1, 1}/2
243: [0, 0, r, c]
244: 	(+0.0000e+00,+0.0000e+00) M	
245: [0, 1, r, c]
246: 	(+7.0711e-01,+0.0000e+00) N	
247: End: {1, 2, 1, 1}/2
248: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{-1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
249: ----------------------------------------
250: Tensor block of rank 4 and 1 symmetries. isZero = +0.
251: Reduced tensor is:
252: Start: {1, 2, 1, 1}/2
253: [0, 0, r, c]
254: 	(+0.0000e+00,+0.0000e+00) O	
255: [0, 1, r, c]
256: 	(+7.0711e-01,+0.0000e+00) P	
257: End: {1, 2, 1, 1}/2
258: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
259: ----------------------------------------
260: Tensor block of rank 4 and 1 symmetries. isZero = +0.
261: Reduced tensor is:
262: Start: {3, 1, 1, 1}/3
263: [0, 0, r, c]
264: 	(+7.0711e-01,+0.0000e+00) Q	
265: [1, 0, r, c]
266: 	(+0.0000e+00,+0.0000e+00) R	
267: [2, 0, r, c]
268: 	(+0.0000e+00,+0.0000e+00) S	
269: End: {3, 1, 1, 1}/3
270: Sym 1 (no  CGC):		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
271: ----------------------------------------
272: Tensor block of rank 4 and 1 symmetries. isZero = +0.
273: Reduced tensor is:
274: Start: {3, 1, 1, 1}/3
275: [0, 0, r, c]
276: 	(+7.0711e-01,+0.0000e+00) T	
277: [1, 0, r, c]
278: 	(+0.0000e+00,+0.0000e+00) U	
279: [2, 0, r, c]
280: 	(+0.0000e+00,+0.0000e+00) V	
281: End: {3, 1, 1, 1}/3
282: Sym 1 (no  CGC):		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
283: ========================================

Still a bit more annoying but still manageable. I've inserted numbers above and added them below to the corresponding entries. Not-stored zero entries are marked with a .. All in all, of 5×4×2×2 = 80 entries, we only store 22 explicitly, of which only 10 are non-zero. The order of the quantum number sectors is somewhat arbitrary, but consistent with the representation on the previous tensor.

284:       <================== Sz=0 ==================> <== Sz=-1 ===> <== Sz=+1 ===>
285: 
286:   Λ   [ -0.5A  .   ] [   0C   .   ] [  0E   .    ] [  .      .  ] [  .    0O   ] 
287:   |   [   .   0.5G ] [   .    0I  ] [  .    0K   ] [  0M     .  ] [  .    .    ] 
288: Sz=0
289:   |   [   0B   .   ] [ -0.5D  .   ] [  1F   .    ] [  .      .  ] [  .    0.7P ] 
290:   V   [   .    0H  ] [   .   0.5J ] [  .    1L   ] [ 0.7M    .  ] [  .    .    ] 
291: 
292: Sz=-1 [   .   0.7T ] [   .   0U   ] [  .    0V   ]    .      .       .    .      
293:       [   .    .   ] [   .    .   ] [  .    .    ]    .      .       .    .
294: 
295: Sz=+1 [   .    .   ] [   .    .   ] [  .    .    ]    .      .       .    .
296:       [  0.7Q  .   ] [   0R   .   ] [  0S   .    ]    .      .       .    .

This tensor can easily be recognised as

297: [ Z 0 0 0 0 ]
298: [ 0 Z 1 + - ]  ± = 1/sqrt(2) S^±
299: [ + 0 0 0 0 ]
300: [ - 0 0 0 0 ]

Continuing, we could sort the other tensor components similarly, but don't.

301: Site: 2
302: ========================================
303: Tensor of rank 4 with the following bases:
304: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <2, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
305: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <3, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
306: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
307: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
308: Composed of the following 6 blocks:
309: ----------------------------------------
310: Tensor block of rank 4 and 1 symmetries. isZero = +0.
311: Reduced tensor is:
312: Start: {2, 3, 1, 1}/6
313: [0, 0, r, c]
314: 	(+1.0000e+00,+0.0000e+00)	
315: [0, 1, r, c]
316: 	(-5.0000e-01,+0.0000e+00)	
317: [0, 2, r, c]
318: 	(+0.0000e+00,+0.0000e+00)	
319: [1, 0, r, c]
320: 	(+0.0000e+00,+0.0000e+00)	
321: [1, 1, r, c]
322: 	(+0.0000e+00,+0.0000e+00)	
323: [1, 2, r, c]
324: 	(-5.0000e-01,+0.0000e+00)	
325: End: {2, 3, 1, 1}/6
326: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
327: ----------------------------------------
328: Tensor block of rank 4 and 1 symmetries. isZero = +0.
329: Reduced tensor is:
330: Start: {2, 3, 1, 1}/6
331: [0, 0, r, c]
332: 	(+1.0000e+00,+0.0000e+00)	
333: [0, 1, r, c]
334: 	(+5.0000e-01,+0.0000e+00)	
335: [0, 2, r, c]
336: 	(+0.0000e+00,+0.0000e+00)	
337: [1, 0, r, c]
338: 	(+0.0000e+00,+0.0000e+00)	
339: [1, 1, r, c]
340: 	(+0.0000e+00,+0.0000e+00)	
341: [1, 2, r, c]
342: 	(+5.0000e-01,+0.0000e+00)	
343: End: {2, 3, 1, 1}/6
344: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
345: ----------------------------------------
346: Tensor block of rank 4 and 1 symmetries. isZero = +0.
347: Reduced tensor is:
348: Start: {1, 3, 1, 1}/3
349: [0, 0, r, c]
350: 	(+0.0000e+00,+0.0000e+00)	
351: [0, 1, r, c]
352: 	(+0.0000e+00,+0.0000e+00)	
353: [0, 2, r, c]
354: 	(+7.0711e-01,+0.0000e+00)	
355: End: {1, 3, 1, 1}/3
356: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{-1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
357: ----------------------------------------
358: Tensor block of rank 4 and 1 symmetries. isZero = +0.
359: Reduced tensor is:
360: Start: {1, 3, 1, 1}/3
361: [0, 0, r, c]
362: 	(+0.0000e+00,+0.0000e+00)	
363: [0, 1, r, c]
364: 	(+0.0000e+00,+0.0000e+00)	
365: [0, 2, r, c]
366: 	(+7.0711e-01,+0.0000e+00)	
367: End: {1, 3, 1, 1}/3
368: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+1}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
369: ----------------------------------------
370: Tensor block of rank 4 and 1 symmetries. isZero = +0.
371: Reduced tensor is:
372: Start: {2, 1, 1, 1}/2
373: [0, 0, r, c]
374: 	(+7.0711e-01,+0.0000e+00)	
375: [1, 0, r, c]
376: 	(+0.0000e+00,+0.0000e+00)	
377: End: {2, 1, 1, 1}/2
378: Sym 1 (no  CGC):		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
379: ----------------------------------------
380: Tensor block of rank 4 and 1 symmetries. isZero = +0.
381: Reduced tensor is:
382: Start: {2, 1, 1, 1}/2
383: [0, 0, r, c]
384: 	(+7.0711e-01,+0.0000e+00)	
385: [1, 0, r, c]
386: 	(+0.0000e+00,+0.0000e+00)	
387: End: {2, 1, 1, 1}/2
388: Sym 1 (no  CGC):		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
389: ========================================
390: 
391: Site: 3
392: ========================================
393: Tensor of rank 4 with the following bases:
394: Leg 1:	Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
395: Leg 2:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>, <2, [U(1)_z:{+0}@1]>, <1, [U(1)_z:{+1}@1]>]
396: Leg 3:	Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
397: Leg 4:	Outgoing Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]
398: Composed of the following 4 blocks:
399: ----------------------------------------
400: Tensor block of rank 4 and 1 symmetries. isZero = +0.
401: Reduced tensor is:
402: Start: {1, 2, 1, 1}/2
403: [0, 0, r, c]
404: 	(+1.0000e+00,+0.0000e+00)	
405: [0, 1, r, c]
406: 	(-5.0000e-01,+0.0000e+00)	
407: End: {1, 2, 1, 1}/2
408: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{-0.5}@1; 	Idx 4 transforms as: U(1)_z:{-0.5}@1
409: ----------------------------------------
410: Tensor block of rank 4 and 1 symmetries. isZero = +0.
411: Reduced tensor is:
412: Start: {1, 2, 1, 1}/2
413: [0, 0, r, c]
414: 	(+1.0000e+00,+0.0000e+00)	
415: [0, 1, r, c]
416: 	(+5.0000e-01,+0.0000e+00)	
417: End: {1, 2, 1, 1}/2
418: Sym 1 (no  CGC):		Idx 1 transforms as: U(1)_z:{+0}@1; 	Idx 2 transforms as: U(1)_z:{+0}@1; 	Idx 3 transforms as: U(1)_z:{+0.5}@1; 	Idx 4 transforms as: U(1)_z:{+0.5}@1
419: ----------------------------------------
420: Tensor block of rank 4 and 1 symmetries. isZero = +0.
421: Reduced tensor is: (+7.0711e-01,+0.0000e+00)
422: Sym 1 (no  CGC):		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
423: ----------------------------------------
424: Tensor block of rank 4 and 1 symmetries. isZero = +0.
425: Reduced tensor is: (+7.0711e-01,+0.0000e+00)
426: Sym 1 (no  CGC):		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
427: ========================================

Working on a very small lattice is fine and obviously useful initially, but those systems can still be diagonalised exactly. The power of tensor product approaches lies in their scaling with lattice size, which is usually polynomial. For states, this scaling is usually \(O(L)\) or \(O(L^2)\) in memory and \(O(L)\) or \(O(L^6)\) in CPU time for most operations for gapped (\(S \sim \mathrm{const}\)) and critical systems (\(S \sim \mathrm{log}(L)\) in one dimension) respectively. The bond dimensions of most operators also saturate with increasing system size or at most grow polynomially, again resulting in approximate \(O(L)\) scaling of memory with system length.

As an example, we can very easily create a longer lattice of length \(L = 100\) and print the bond dimension of Hp Hz +:

428: $ syten-sql-mps-spin --sym u1 -l 100 -z 0.5 -o sp100.lat
429: I: Creating a 100×1 spin lattice in file sp100.lat, using mapping variant Z and shift 0 with u1 symmetries.
430:   Ordering is defined by new[old]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99].
431: ....v....x....v....x....v....x....v....x....v....l....v....x....v....x....v....x....v....x....v....c
432: 2017.08.09T16:17:46.672 1.82807e+00 I Working on operator Hp
433: 2017.08.09T16:17:47.125 2.28150e+00 I Working on operator Hp_xp
434: 2017.08.09T16:17:47.126 2.28163e+00 I Working on operator Hz
435: 2017.08.09T16:17:47.732 2.88010e+00 I Working on operator Hz_s
436: 2017.08.09T16:17:48.264 3.41226e+00 I Working on operator Hz_xp
437: Saving lattice…
438: Saved lattice into file sp100.lat.
439: $ syten-info -qb sp100.lat:'Hz Hp +'
440:      (0,1)         3         4         4         1         2
441:      (1,2)         3         5         5         2         3
442:      (2,3)         3         5         5         2         3
443:      (3,4)         3         5         5         2         3
444: […snip…]
445:    (95,96)         3         5         5         2         3
446:    (96,97)         3         5         5         2         3
447:    (97,98)         3         5         5         2         3
448:    (98,99)         3         4         4         1         2

Here, we can clearly see that the bond dimension of the operator Hp Hz + is independent of the length of the spin chain.

2.3 Generating a random initial state

Use syten-random with the lattice in question:

449: $ syten-random -l sp4.lat -o rnd
450: Seed used: 4270940170
451: Loaded MP lattice file.
452: Parsing sector…
453: Creating a sampled state transforming as Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>] with m = 10 and m_b = 10 on 4 sites.
454: ....v....x
455: Adding generated states…
456: Truncating states…
457: Saving state to rnd

By default, a state with quantum number 0 will be created by applying random operators in random places until the overall state transforms as desired:

458: $ syten-info rnd
459: 2017.08.09T16:30:23.809 2.09170e-02 I Loaded MP state
460: 2017.08.09T16:30:23.809 2.10140e-02 I MP State on lattice of 4 sites.
461: 2017.08.09T16:30:23.809 2.10670e-02 I Vacuum basis (right edge): [<1, [U(1)_z:{+0}@1]>]
462: 2017.08.09T16:30:23.809 2.11110e-02 I Sector basis (left edge):  [<1, [U(1)_z:{+0}@1]>]
463: 2017.08.09T16:30:23.809 2.11320e-02 I Physical basis on site 1:  [<1, [U(1)_z:{+0.5}@1]>, <1, [U(1)_z:{-0.5}@1]>]

Due to the small lattice size, its bond dimension will be small (but also maximal, d(L/2)):

464: $ syten-info rnd -qb
465: (0,1)         2         2         2         1         1
466: (1,2)         3         4         4         1         2
467: (2,3)         2         2         2         1         1

Possible to specify a different quantum number sector:

468: $ syten-random -l sp4.lat -o rnd -s 2
469: Seed used: 3781107722
470: Loaded MP lattice file.
471: Parsing sector…
472: Creating a sampled state transforming as Outgoing Basis (1 Syms): [<1, [U(1)_z:{+2}@1]>] with m = 10 and m_b = 10 on 4 sites.
473: ....v....x
474: Adding generated states…
475: Truncating states…
476: Saving state to rnd

No checking is done, if sector is unreachable, this will not terminate:

477: $ syten-random -l sp4.lat -o rnd -s 2.5
478: Seed used: 1320226826
479: Loaded MP lattice file.
480: Parsing sector…
481: Creating a sampled state transforming as Outgoing Basis (1 Syms): [<1, [U(1)_z:{+2.5}@1]>] with m = 10 and m_b = 10 on 4 sites.
482: ^C

2.4 Expectation values

We can already calculate expectation values of, e.g., Hz or Hp:

483: $ syten-expectation rnd sp4.lat:Hp -q
484: (0,0)
485: $ syten-expectation rnd sp4.lat:Hz -q
486: (0.75,3.24651169596159e-19)
487: $ syten-expectation rnd sp4.lat:'Hz Hp +' -q
488: (0.75,8.11462711676569e-20)
489: $ syten-localexpectation rnd sp4.lat sz
490: 2017.08.09T16:36:42.963 2.93400e-02 I Loading first state…
491: 2017.08.09T16:36:42.964 3.03560e-02 I Loading lattice…
492: 2017.08.09T16:36:42.975 3.75920e-02 N Overlap: (1,0)
493: Site                           Real part            Imaginary part
494: 0                 +5.000000000000000e-01    +0.000000000000000e+00 
495: 1                 +5.000000000000000e-01    +0.000000000000000e+00 
496: 2                 +5.000000000000000e-01    +0.000000000000000e+00 
497: 3                 +5.000000000000000e-01    +0.000000000000000e+00

Q: Why is this all +½?

2.5 Operator application

498: $ syten-apply-op -i rnd -l sp4.lat:'sm 0 @' -o rnd-sm0
499: 2017.08.09T16:37:42.179 3.20640e-02 I Loaded MP state file 'rnd'.
500: 2017.08.09T16:37:42.179 3.21770e-02 I Loading and parsing MP operator 'sp4.lat:sm 0 @'…
501: 2017.08.09T16:37:42.182 3.53640e-02 I Applying operator naively…
502: 2017.08.09T16:37:42.187 3.94360e-02 I Truncating resulting state to δs = 2.22045e-15 and m = 4294967295…
503: 2017.08.09T16:37:42.188 4.04680e-02 N Incurred error: 0.
504: 2017.08.09T16:37:42.188 4.05190e-02 I Saving to file 'rnd-sm0'…
505: $ syten-localexpectation rnd-sm0 sp4.lat sz
506: 2017.08.09T16:37:47.369 1.95250e-02 I Loading first state…
507: 2017.08.09T16:37:47.370 2.07600e-02 I Loading lattice…
508: 2017.08.09T16:37:47.380 2.78030e-02 N Overlap: (1,0)
509: Site                           Real part            Imaginary part
510: 0                 -5.000000000000000e-01    +0.000000000000000e+00 
511: 1                 +5.000000000000000e-01    +0.000000000000000e+00 
512: 2                 +5.000000000000000e-01    +0.000000000000000e+00 
513: 3                 +5.000000000000000e-01    +0.000000000000000e+00
514: $ syten-apply-op rnd sp4.lat:Hz -o rnd-Hz
515: 2017.08.09T16:38:40.246 2.02070e-02 I Loaded MP state file 'rnd'.
516: 2017.08.09T16:38:40.246 2.02910e-02 I Loading and parsing MP operator 'sp4.lat:Hz'…
517: 2017.08.09T16:38:40.248 2.24250e-02 I Applying operator naively…
518: 2017.08.09T16:38:40.254 2.76420e-02 I Truncating resulting state to δs = 2.22045e-15 and m = 4294967295…
519: 2017.08.09T16:38:40.254 2.85600e-02 N Incurred error: 0.
520: 2017.08.09T16:38:40.255 2.86010e-02 I Saving to file 'rnd-Hz'…
521: $ syten-norm rnd
522: Loaded MP state, calculating norm…
523: 1
524: $ syten-norm rnd-Hz
525: Loaded MP state, calculating norm…
526: 0.75
527: $ syten-overlap rnd-Hz rnd
528: 2017.08.09T16:39:27.204 3.25470e-02 N Loaded first MP state.
529: 2017.08.09T16:39:27.204 3.28000e-02 N Loaded second MP state.
530: 2017.08.09T16:39:27.204 3.28260e-02 N Calculating overlap…
531: (0.7500000000000002,0)

The maximal bond dimension will grow with every operator application, unless the state is an eigenstate of the operator or the operator has bond dimension 1:

532: $ syten-random -l sp100.lat -o rnd-100 -s 0
533: Seed used: 1661235210
534: Loaded MP lattice file.
535: Parsing sector…
536: Creating a sampled state transforming as Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>] with m = 10 and m_b = 10 on 100 sites.
537: ....v....x
538: Adding generated states…
539: Truncating states…
540: Saving state to rnd-100
541: $ syten-info --get-max-total -q rnd-100
542: 10
543: $ syten-apply-op -i rnd-100 -o rnd-100-H -l sp100.lat:'Hz Hp +'
544: 2017.08.09T16:43:31.456 3.18010e-02 I Loaded MP state file 'rnd-100'.
545: 2017.08.09T16:43:31.456 3.18940e-02 I Loading and parsing MP operator 'sp100.lat:Hz Hp +'…
546: 2017.08.09T16:43:31.594 1.68797e-01 I Applying operator naively…
547: 2017.08.09T16:43:31.641 2.15027e-01 I Truncating resulting state to δs = 2.22045e-15 and m = 4294967295…
548: 2017.08.09T16:43:31.702 2.76756e-01 N Incurred error: 0.
549: 2017.08.09T16:43:31.702 2.76809e-01 I Saving to file 'rnd-100-H'…
550: $ syten-info --get-max-total -q rnd-100-H
551: 30

3 MPS-DMRG application to a spin chain: the power of symmetries

In the previous section, we have used a \(\mathrm{U}(1)_{S^z}\)-symmetric spin chain. Symmetries allow the labelling of states according to their quantum number and introduce additional sparsity into the MPS and MPO tensors. Using them can greatly speed up a computation. As an example, generate three lattices with no symmetry, \(\mathrm{U}(1)_{S^z}\) and \(\mathrm{SU}(2)_{\mathrm{Spin}}\) symmetry respectively:

552: $ syten-sql-mps-spin -l 30 -o sp-nil --sym nil -q
553: $ syten-sql-mps-spin -l 30 -o sp-u1 --sym u1 -q
554: $ syten-sql-mps-spin -l 30 -o sp-su2 --sym su2 -q

Then, generate three random states:

555: $ for sym in nil u1 su2; do syten-random -l sp-${sym} -o rnd-${sym} -q; done

And prepare by running DMRG on the "no-symmetries" state. DMRG requires you to specify a Hamiltonian, an initial state, optionally an output prefix and a staging configuration which specifies how many sweeps to do with which bond dimension, truncation error etc. For now, we will do up to twenty sweeps at m=200 without additional truncation. Take note that the ‘nil’-lattice does not define a global H or Hp and Hz but only H{x,y,z}, which we have to add together:

556: $ syten-dmrg -l sp-nil:'Hx Hy Hz + +' -i rnd-nil -s "(m 200 x 20 t 0)" -o dmrg-nil
557: […snip…]
558: 2017.08.09T17:30:38.959 6.40253e+01    1  ←    1    7     1  -1.31113557586031e+01   1      nan  -1.31113557586031e+01  0.00e+00  1.02e+00     1     1     6       6     6    invSubsp[n=2.3e-10]
559: […snip…]
560: 2017.08.09T17:30:38.959 6.40254e+01 N  1 Completed right-to-left (←) sweep 7 of stage 1.
561: 2017.08.09T17:30:38.959 6.40254e+01 N  1 Convergence: Converged: Energy gain threshold reached: -6.63863e-15 < 2.22045e-15
562: 2017.08.09T17:30:38.959 6.40254e+01 N  1 Report: (Converged) Energy gain threshold reached at -6.63863e-15 < 2.22045e-15
563: 2017.08.09T17:30:38.959 6.40254e+01 N  M All workers converged.
564: 2017.08.09T17:30:38.959 6.40255e+01 N  1 Extrapolation info (minimal energy, maximal truncation):
565: 2017.08.09T17:30:38.959 6.40255e+01 N  1 Energy/Trunc: -13.111355758603196 0
566: 2017.08.09T17:30:38.959 6.40255e+01 N  M Completed stage 1…
567: 2017.08.09T17:30:38.959 6.40255e+01 N  M Saving state…
568: 2017.08.09T17:30:38.959 6.40255e+01 N  M Copying state & right-normalising…
569: 2017.08.09T17:30:39.396 6.44630e+01 N  M Saved state to file dmrg-nil_1_7.state

After 64 seconds and after the seventh sweep, the calculation has completed. We can compare this to the \(\mathrm{U}(1)_{S^z}\)-symmetric case:

570: $ syten-dmrg -l sp-u1:'Hp Hz +' -i rnd-u1 -s "(m 200 x 20 t 0)" -o dmrg-u1
571: […snip…]
572: 2017.08.09T17:32:05.218 1.24163e+01    1  →    1   12    28  -1.31113557585980e+01   1      nan  -1.31113557585980e+01  0.00e+00  8.22e+01     6     4     8       8     3    invSubsp[n=4.5e-11]
573: […snip…]
574: 2017.08.09T17:32:05.218 1.24164e+01 N  1 Completed left-to-right (→) sweep 12 of stage 1.
575: 2017.08.09T17:32:05.218 1.24164e+01 N  1 Convergence: Converged: Energy gain threshold reached: -4.06447e-14 < 2.22045e-15
576: 2017.08.09T17:32:05.218 1.24164e+01 N  1 Report: (Converged) Energy gain threshold reached at -4.06447e-14 < 2.22045e-15
577: 2017.08.09T17:32:05.218 1.24164e+01 N  M All workers converged.
578: 2017.08.09T17:32:05.218 1.24164e+01 N  1 Extrapolation info (minimal energy, maximal truncation):
579: 2017.08.09T17:32:05.218 1.24164e+01 N  1 Energy/Trunc: -13.111355758601112 2.5660028091299701e-10
580: 2017.08.09T17:32:05.218 1.24164e+01 N  M Completed stage 1…
581: 2017.08.09T17:32:05.218 1.24165e+01 N  M Saving state…
582: 2017.08.09T17:32:05.218 1.24165e+01 N  M Copying state & right-normalising…
583: 2017.08.09T17:32:05.272 1.24699e+01 N  M Saved state to file dmrg-u1_1_12.state

Due to a different initial state, the calculation now took 12 sweeps, but only 12 seconds to complete them! Going on to the full \(\mathrm{SU}(2)_{\mathrm{Spin}}\)-symmetric case:

584: $ syten-dmrg -l sp-su2:'H' -i rnd-su2 -s "(m 200 x 20 t 0)" -o dmrg-su2
585: […snip…]
586: 2017.08.09T17:33:37.329 7.78297e+00    1  ←    1    7     1  -1.31113557586031e+01   1      nan  -1.31113557586031e+01  0.00e+00  9.18e-01     3     2     4      10     3    invSubsp[n=3.3e-09]
587: […snip…]
588: 2017.08.09T17:33:37.329 7.78308e+00 N  1 Completed right-to-left (←) sweep 7 of stage 1.
589: 2017.08.09T17:33:37.329 7.78311e+00 N  1 Convergence: Converged: Energy gain threshold reached: -4.06447e-16 < 2.22045e-15
590: 2017.08.09T17:33:37.329 7.78312e+00 N  1 Report: (Converged) Energy gain threshold reached at -4.06447e-16 < 2.22045e-15
591: 2017.08.09T17:33:37.329 7.78314e+00 N  M All workers converged.
592: 2017.08.09T17:33:37.329 7.78315e+00 N  1 Extrapolation info (minimal energy, maximal truncation):
593: 2017.08.09T17:33:37.329 7.78316e+00 N  1 Energy/Trunc: -13.111355758603127 0
594: 2017.08.09T17:33:37.329 7.78318e+00 N  M Completed stage 1…
595: 2017.08.09T17:33:37.329 7.78319e+00 N  M Saving state…
596: 2017.08.09T17:33:37.329 7.78320e+00 N  M Copying state & right-normalising…
597: 2017.08.09T17:33:37.429 7.88359e+00 N  M Saved state to file dmrg-su2_1_7.state

Just 7.8 seconds now! When checking the output carefully, you will also notice that there was also a major change in the last columns (nil, \(\mathrm{U}(1)_{S^z}\) and \(\mathrm{SU}(2)_{\mathrm{Spin}}\) atop of each other):

598: $ grep '→    1    6    15' *i1log | column -t
599:                                                                                                                       size of the largest dense block (O(m³) scaling) +
600:                                                                                   bond dimension 'seen' by the phsical system due to inner multiplicity of SU(2) +    |
601:                                                                                                     effective bond dimension, the value specified using 'm' +    |    |
602:                                                                                                                                                             v    v    v
603: dmrg-nil.i1log:2017.08.09T17:30:21.136  4.62016e+01  1  →  1  6  15  -1.31113557586031e+01  8  0.0e+00  -1.31113557586031e+01  0.00e+00  1.02e+00  1   1  200  200  200  maxIter
604: dmrg-su2.i1log:2017.08.09T17:33:35.861  6.31513e+00  1  →  1  6  15  -1.31113557586031e+01  1  nan      -1.31113557586031e+01  0.00e+00  9.14e-01  10  6  200  748  76   invSubsp[n=7.6e-09]
605: dmrg-u1.i1log:2017.08.09T17:31:56.821   4.01858e+00  1  →  1  6  15  -1.31113557586030e+01  8  5.8e-04  -1.31113557586030e+01  1.71e-05  1.00e+02  11  6  200  200  71   maxIter

Implementing the inner multiplicities from \(\mathrm{SU}(2)_{\mathrm{Spin}}\) allows representing more states at the same effective bond dimension m. Implementing any symmetry allows splitting the dense blocks into sub-blocks of much smaller size. Compare also the file size of each generated wavefunction:

606: $ ls -lh *.state
607: -rw-rw-rw- 1 C.Hubig ls-schollwoeck  22M 2017-08-09 17:30 dmrg-nil_1_7.state
608: -rw-rw-rw- 1 C.Hubig ls-schollwoeck 2.9M 2017-08-09 17:33 dmrg-su2_1_7.state
609: -rw-rw-rw- 1 C.Hubig ls-schollwoeck 4.5M 2017-08-09 17:32 dmrg-u1_1_12.state

We can also check that we have found an eigenstate of the Hamiltonian by calculating the variance, which should be small:

610: $ syten-expectation -qrv dmrg-su2_1_7.state sp-su2:H
611: ⟨O⟩:          -13.1113557586031
612: ⟨O²⟩:         171.907649828656
613: |⟨O²⟩-⟨O⟩²|:  2.8421709430404e-14

3.1 Efficient calculation of correlators \(\langle \hat S^z_i \hat S^z_{15} \rangle\)

First, apply Sz15 to the ground state:

614: $ syten-apply-op -i dmrg-nil_1_7.state -o gs-nil-sz15 -l sp-nil:'sz 15 @'

Then, calculate the local expectation values \(\langle \hat S^z_{15} | \hat S^z_i |0\rangle\), taking care not to normalise by \(\langle \hat S^z_{15}|0\rangle\) (via -s):

615: $ syten-localexpectation gs-nil-sz15 sp-nil sz dmrg-nil_1_7.state -s
616: 2017.08.09T17:42:30.137 1.94260e-02 I Loading first state…
617: 2017.08.09T17:42:30.255 1.37436e-01 I Loading second state…
618: 2017.08.09T17:42:30.340 2.22686e-01 I Loading lattice…
619: 2017.08.09T17:42:31.288 1.16372e+00 N Overlap: assumed 1
620: Site                           Real part            Imaginary part
621: 0                 -1.179568834873305e-02    +9.107298248878237e-18 
622: 1                 +7.697983540194697e-03    -4.662069341687669e-18 
623: 2                 -1.477432218289132e-02    +1.127570259384925e-17 
624: 3                 +1.104464016318039e-02    -9.974659986866641e-18 
625: 4                 -1.784312419397555e-02    -6.938893903907228e-18 
626: 5                 +1.425029623205436e-02    -7.806255641895632e-18 
627: 6                 -2.174962728844476e-02    +1.821459649775647e-17 
628: 7                 +1.820037388176275e-02    -5.637851296924623e-18 
629: 8                 -2.750091819303292e-02    -1.734723475976807e-18 
630: 9                 +2.401872368495444e-02    +2.168404344971009e-17 
631: 10                -3.743469033826867e-02    -1.387778780781446e-17 
632: 11                +3.445571496313890e-02    -3.035766082959412e-18 
633: 12                -5.983807317348300e-02    +8.673617379884035e-18 
634: 13                +6.062733073252365e-02    -2.775557561562891e-17 
635: 14                -1.686782160834168e-01    -8.673617379884035e-17 
636: 15                +2.500000000000006e-01    +2.411807732694005e-16 
637: 16                -1.268241155518576e-01    -1.457167719820518e-16 
638: 17                +6.074810072697510e-02    +3.339342691255354e-17 
639: 18                -4.068663796789525e-02    -3.122502256758253e-17 
640: 19                +3.476377533148065e-02    +2.005774019098183e-17 
641: 20                -2.421376298714984e-02    +2.862293735361732e-17 
642: 21                +2.455178413838127e-02    +1.501620008892424e-17 
643: 22                -1.704824587427810e-02    +1.435212625827686e-17 
644: 23                +1.900562396660125e-02    -2.982826523756019e-17 
645: 24                -1.280714303281470e-02    +2.452075679006474e-17 
646: 25                +1.539965929320117e-02    -2.753873518113181e-17 
647: 26                -9.677886949512262e-03    +2.145216818039922e-17 
648: 27                +1.266461730534876e-02    -2.825701912040346e-17 
649: 28                -6.641172955762460e-03    +1.500059615384904e-17 
650: 29                +1.008500116172125e-02    -2.404908346213956e-17

You can see that there are, even at L=30, still finite-size effects. The SU(2)-symmetric case has the same problem. Note that we can only calculate \(\langle \hat S_{i} \cdot \hat S_{j}\rangle =\) \(\langle \hat S^x_i \cdot \hat S^x_j\rangle + \langle \hat S^y_i \cdot \hat S^y_j \rangle + \langle \hat S^z_i \cdot \hat S^z_j\rangle\), not individual components, so we divide the real part by three with awk to ease comparison with the calculation above:

651: $ syten-apply-op -i dmrg-su2_1_7.state -o gs-su2-s15 -l sp-su2:'s 15 @'
652: $ syten-localexpectation gs-su2-s15 sp-su2 s dmrg-su2_1_7.state -s -q | awk '{print $0, $2/3.}' | column -t
653: 0   -3.538706291169211e-02  +0.000000000000000e+00  -0.0117957
654: 1   +2.309394917897722e-02  +0.000000000000000e+00  0.00769798
655: 2   -4.432296487907535e-02  +0.000000000000000e+00  -0.0147743
656: 3   +3.313391916665966e-02  +0.000000000000000e+00  0.0110446
657: 4   -5.352937141815489e-02  +0.000000000000000e+00  -0.0178431
658: 5   +4.275088764130345e-02  +0.000000000000000e+00  0.0142503
659: 6   -6.524888124362017e-02  +0.000000000000000e+00  -0.0217496
660: 7   +5.460112091136803e-02  +0.000000000000000e+00  0.0182004
661: 8   -8.250275445571692e-02  +0.000000000000000e+00  -0.0275009
662: 9   +7.205617057525550e-02  +0.000000000000000e+00  0.0240187
663: 10  -1.123040713172472e-01  +0.000000000000000e+00  -0.0374347
664: 11  +1.033671446130366e-01  +0.000000000000000e+00  0.0344557
665: 12  -1.795142202935994e-01  +0.000000000000000e+00  -0.0598381
666: 13  +1.818819920720938e-01  +0.000000000000000e+00  0.0606273
667: 14  -5.060346501414829e-01  +0.000000000000000e+00  -0.168678
668: 15  +7.500000000000026e-01  +0.000000000000000e+00  0.25
669: 16  -3.804723448403237e-01  +0.000000000000000e+00  -0.126824
670: 17  +1.822443020713638e-01  +0.000000000000000e+00  0.0607481
671: 18  -1.220599131020208e-01  +0.000000000000000e+00  -0.0406866
672: 19  +1.042913258088060e-01  +0.000000000000000e+00  0.0347638
673: 20  -7.264128841647158e-02  +0.000000000000000e+00  -0.0242138
674: 21  +7.365535217206616e-02  +0.000000000000000e+00  0.0245518
675: 22  -5.114473720429411e-02  +0.000000000000000e+00  -0.0170482
676: 23  +5.701687162013721e-02  +0.000000000000000e+00  0.0190056
677: 24  -3.842142875953879e-02  +0.000000000000000e+00  -0.0128071
678: 25  +4.619897756710443e-02  +0.000000000000000e+00  0.0153997
679: 26  -2.903366057313903e-02  +0.000000000000000e+00  -0.00967789
680: 27  +3.799385161981647e-02  +0.000000000000000e+00  0.0126646
681: 28  -1.992351869179369e-02  +0.000000000000000e+00  -0.00664117
682: 29  +3.025500323018018e-02  +0.000000000000000e+00  0.010085

4 MPS-DMRG in unconnected sub spaces [extra]

When applying DMRG to matrix-product states, one has to keep in mind that the method boils down to a very fancy version of the Lanczos method. This is also true if one of the typical perturbation terms such as subspace expansion or the density matrix perturbation are used. As such, if for the initial state \(|\psi\rangle\) and Hamiltonian \(\hat H\), there is no \(n\) for which \(\hat H^n |\psi\rangle\) has non-zero overlap with the ground state, the method is doomed to fail. In typical model applications, this case is extremely rare. In applications to real materials or embedding problems, however, it is not untypical to encounter cases where this occurs – the typical example is a Hamiltonian which does not allow for particles to hop between parts of the lattice and hence cannot re-distribute them there. This can be extended to allow for pair-hopping but not exchange of a single particle, leading to preservation of sub-lattice parity. Additionally, the "unconnected" parts don't even have to be adjacent, they might as well be spread throughout the system.

One ‘model’ example where this situation does occur is the Fermi-Hubbard model with just a local interaction \(U\) but no hopping between sites (\(t=0\)):

$$\hat H = U \hat H_U + t \hat H_t = U \sum_{i} \hat n_i^2 - n_i + t \sum_i \hat c^\dagger_i \cdot c_{i+1} + \mathrm{h.c.}$$

We can easily generate a chain of twenty sites and a random initial state with 20 electrons (20,0 stands for twenty electrons (\(N=1\)) with total spin \(S=0\)):

683: $ syten-sql-mps-fermi-hubbard -l 20 -w 1 -o fh
684: $ syten-random -l fh -s 20,0 -o rnd

Now running DMRG on this initial state will prove disappointing (we can try with both single-site and two-site DMRG):

685: $ syten-dmrg -l fh:Hd -i rnd -s "(m 100 v DMRG3S)" -o dmrg3s -f dmrg3s
686: $ syten-dmrg -l fh:Hd -i rnd -s "(m 100 v 2DMRG )" -o 2dmrg -f 2dmrg
687: $ paste <(syten-localexpectation -qr 2dmrg fh n) <(syten-localexpectation -qr dmrg3s fh n | awk '{print $2}') | column -t
688: 0   +1.000000000000000e+00  +1.000000000000000e+00
689: 1   +1.000000000000000e+00  +1.000000000000000e+00
690: 2   +5.000000000000057e-01  +5.000000000000063e-01
691: 3   +5.000000000000057e-01  +5.000000000000063e-01
692: 4   +1.500000000000006e+00  +1.500000000000006e+00
693: 5   +5.000000000000057e-01  +5.000000000000063e-01
694: 6   +4.999999999999944e-01  +4.999999999999937e-01
695: 7   +1.499999999999994e+00  +1.499999999999994e+00
696: 8   +1.000000000000000e+00  +1.000000000000000e+00
697: 9   +2.000000000000000e+00  +2.000000000000000e+00
698: 10  +1.500000000000006e+00  +1.500000000000006e+00
699: 11  +4.999999999999944e-01  +4.999999999999937e-01
700: 12  +1.000000000000000e+00  +1.000000000000000e+00
701: 13  +1.000000000000000e+00  +1.000000000000000e+00
702: 14  +1.500000000000006e+00  +1.500000000000006e+00
703: 15  +5.000000000000057e-01  +5.000000000000063e-01
704: 16  +9.999999999999888e-01  +9.999999999999875e-01
705: 17  +9.999999999999888e-01  +9.999999999999875e-01
706: 18  +1.000000000000011e+00  +1.000000000000013e+00
707: 19  +9.999999999999888e-01  +9.999999999999875e-01

Evidently, DMRG failed to re-distribute the particles through the system. To allow for this, we will add a small hopping term Ht, which will be multiplied by a factor of 0.5 at every sweep and hence slowly taken to zero. This should allow for the particles to ‘settle’ in their optimal positions.

708: $ syten-dmrg -l fh:Hd -l fh:Ht -i rnd -s "(m 100 v DMRG3S cs 1 0.5 l 10 ) (cs 1 0 x 4)" -o pert -f pert

Note that we pass in two Hamiltonians and then add the rather arcane "cs 1 0.5" to the staging configuration. This says to multiply the 1-st Hamiltonian (zero-indexed, hence fh:Ht) by 0.5 at the end of every sweep. Since the energy will now not necessarily be monotonous (it may go up if the perturbation changes), we specify a minimal number of 10 sweeps to avoid DMRG exiting early. In a second stage, marked with parantheses, we set the coefficient to exactly zero (cs 1 0) and limit the maximal number of sweeps to four (x 4).

The result looks much better:

709: $ syten-localexpectation pert fh n -qr
710: 0                 +1.000000000000000e+00 
711: 1                 +1.000000000000000e+00 
712: 2                 +1.000000000000000e+00 
713: 3                 +1.000000000000000e+00 
714: 4                 +1.000000000000000e+00 
715: 5                 +1.000000000000000e+00 
716: 6                 +1.000000000000000e+00 
717: 7                 +1.000000000000000e+00 
718: 8                 +1.000000000000000e+00 
719: 9                 +1.000000000000000e+00 
720: 10                +1.000000000000000e+00 
721: 11                +1.000000000000000e+00 
722: 12                +1.000000000000000e+00 
723: 13                +1.000000000000000e+00 
724: 14                +1.000000000000000e+00 
725: 15                +1.000000000000000e+00 
726: 16                +1.000000000000000e+00 
727: 17                +1.000000000000000e+00 
728: 18                +1.000000000000000e+00 
729: 19                +1.000000000000000e+00

Additionally, the energy is now much lower and certainly the ground-state energy of \(\hat H_U\):

730: $ for i in 2dmrg dmrg3s pert; do printf "%10s: %10g\n" "${i}" "$(syten-expectation -qr ${i} fh:Hd)"; done
731:  2dmrg:          5
732: dmrg3s:          5
733:   pert:          0

In real-world applications, adding a small hopping term throughout the system which is slowly taken to zero is surprisingly effective at overcoming convergence problems.

5 Time evolution: spectral functions

Time evolution is possible using the matrix-product framework. However, due to the usually quickly growing entanglement in the system, long real-time scales are difficult to achieve or require extraordinary computational resources. Imaginary time evolution is usually less problematic. In the following, we will calculate the real-time correlator \(\langle 0| S^+_i(t) S^-_{32}(0)|0\rangle\) where \(|0\rangle\) is the ground state of a uniform nearest-neighbour spin chain Hamiltonian with an additional, staggered external field of strength 0.3. We will be using a second-order TEBD time evolution method as presented in the lecture. Since this method requires the decomposition of the Hamiltonian into mutually commuting terms, some work was already done preparing an appropriate description of the decomposed Hamiltonian (which is, however, limited to systems of size 65).

First, we prepare the lattice and the random state and run a small DMRG calculation to get the ground state. The staggered field has directions ↑↓↑…↓↑↓↑, so our state should have quantum number -½:

734: $ syten-sql-mps-spin --sym u1 -l 65 -o sp
735: $ syten-random -l sp -o rnd -s -0.5
736: $ syten-dmrg -l sp:'Hp Hz + Hz_s 0.3 * +' -i rnd -s "(m 200 t 0)" -o dmrg -f gs

Then, create the time evolution file containing the operator. Each line of this file names a category for the line, two sites on which the operator acts and a standard MPO description (this file is also kept at /project/summerschool/hubig-hands-on/tebd-swappable-evolution on the tutorial computers):

A 0 1 sz 0 @ sz 1 @ * sp 0 @ sm 1 @ * sm 0 @ sp 1 @ * + 0.5 * + sz 0 @ 0.3 * +
A 2 3 sz 2 @ sz 3 @ * sp 2 @ sm 3 @ * sm 2 @ sp 3 @ * + 0.5 * + sz 2 @ 0.3 * +
A 4 5 sz 4 @ sz 5 @ * sp 4 @ sm 5 @ * sm 4 @ sp 5 @ * + 0.5 * + sz 4 @ 0.3 * +
A 6 7 sz 6 @ sz 7 @ * sp 6 @ sm 7 @ * sm 6 @ sp 7 @ * + 0.5 * + sz 6 @ 0.3 * +
A 8 9 sz 8 @ sz 9 @ * sp 8 @ sm 9 @ * sm 8 @ sp 9 @ * + 0.5 * + sz 8 @ 0.3 * +
A 10 11 sz 10 @ sz 11 @ * sp 10 @ sm 11 @ * sm 10 @ sp 11 @ * + 0.5 * + sz 10 @ 0.3 * +
A 12 13 sz 12 @ sz 13 @ * sp 12 @ sm 13 @ * sm 12 @ sp 13 @ * + 0.5 * + sz 12 @ 0.3 * +
A 14 15 sz 14 @ sz 15 @ * sp 14 @ sm 15 @ * sm 14 @ sp 15 @ * + 0.5 * + sz 14 @ 0.3 * +
A 16 17 sz 16 @ sz 17 @ * sp 16 @ sm 17 @ * sm 16 @ sp 17 @ * + 0.5 * + sz 16 @ 0.3 * +
A 18 19 sz 18 @ sz 19 @ * sp 18 @ sm 19 @ * sm 18 @ sp 19 @ * + 0.5 * + sz 18 @ 0.3 * +
A 20 21 sz 20 @ sz 21 @ * sp 20 @ sm 21 @ * sm 20 @ sp 21 @ * + 0.5 * + sz 20 @ 0.3 * +
A 22 23 sz 22 @ sz 23 @ * sp 22 @ sm 23 @ * sm 22 @ sp 23 @ * + 0.5 * + sz 22 @ 0.3 * +
A 24 25 sz 24 @ sz 25 @ * sp 24 @ sm 25 @ * sm 24 @ sp 25 @ * + 0.5 * + sz 24 @ 0.3 * +
A 26 27 sz 26 @ sz 27 @ * sp 26 @ sm 27 @ * sm 26 @ sp 27 @ * + 0.5 * + sz 26 @ 0.3 * +
A 28 29 sz 28 @ sz 29 @ * sp 28 @ sm 29 @ * sm 28 @ sp 29 @ * + 0.5 * + sz 28 @ 0.3 * +
A 30 31 sz 30 @ sz 31 @ * sp 30 @ sm 31 @ * sm 30 @ sp 31 @ * + 0.5 * + sz 30 @ 0.3 * +
A 32 33 sz 32 @ sz 33 @ * sp 32 @ sm 33 @ * sm 32 @ sp 33 @ * + 0.5 * + sz 32 @ 0.3 * +
A 34 35 sz 34 @ sz 35 @ * sp 34 @ sm 35 @ * sm 34 @ sp 35 @ * + 0.5 * + sz 34 @ 0.3 * +
A 36 37 sz 36 @ sz 37 @ * sp 36 @ sm 37 @ * sm 36 @ sp 37 @ * + 0.5 * + sz 36 @ 0.3 * +
A 38 39 sz 38 @ sz 39 @ * sp 38 @ sm 39 @ * sm 38 @ sp 39 @ * + 0.5 * + sz 38 @ 0.3 * +
A 40 41 sz 40 @ sz 41 @ * sp 40 @ sm 41 @ * sm 40 @ sp 41 @ * + 0.5 * + sz 40 @ 0.3 * +
A 42 43 sz 42 @ sz 43 @ * sp 42 @ sm 43 @ * sm 42 @ sp 43 @ * + 0.5 * + sz 42 @ 0.3 * +
A 44 45 sz 44 @ sz 45 @ * sp 44 @ sm 45 @ * sm 44 @ sp 45 @ * + 0.5 * + sz 44 @ 0.3 * +
A 46 47 sz 46 @ sz 47 @ * sp 46 @ sm 47 @ * sm 46 @ sp 47 @ * + 0.5 * + sz 46 @ 0.3 * +
A 48 49 sz 48 @ sz 49 @ * sp 48 @ sm 49 @ * sm 48 @ sp 49 @ * + 0.5 * + sz 48 @ 0.3 * +
A 50 51 sz 50 @ sz 51 @ * sp 50 @ sm 51 @ * sm 50 @ sp 51 @ * + 0.5 * + sz 50 @ 0.3 * +
A 52 53 sz 52 @ sz 53 @ * sp 52 @ sm 53 @ * sm 52 @ sp 53 @ * + 0.5 * + sz 52 @ 0.3 * +
A 54 55 sz 54 @ sz 55 @ * sp 54 @ sm 55 @ * sm 54 @ sp 55 @ * + 0.5 * + sz 54 @ 0.3 * +
A 56 57 sz 56 @ sz 57 @ * sp 56 @ sm 57 @ * sm 56 @ sp 57 @ * + 0.5 * + sz 56 @ 0.3 * +
A 58 59 sz 58 @ sz 59 @ * sp 58 @ sm 59 @ * sm 58 @ sp 59 @ * + 0.5 * + sz 58 @ 0.3 * +
A 60 61 sz 60 @ sz 61 @ * sp 60 @ sm 61 @ * sm 60 @ sp 61 @ * + 0.5 * + sz 60 @ 0.3 * +
A 62 63 sz 62 @ sz 63 @ * sp 62 @ sm 63 @ * sm 62 @ sp 63 @ * + 0.5 * + sz 62 @ 0.3 * +
B 1 2 sz 1 @ sz 2 @ * sp 1 @ sm 2 @ * sm 1 @ sp 2 @ * + 0.5 * + sz 1 @ -0.3 * +
B 3 4 sz 3 @ sz 4 @ * sp 3 @ sm 4 @ * sm 3 @ sp 4 @ * + 0.5 * + sz 3 @ -0.3 * +
B 5 6 sz 5 @ sz 6 @ * sp 5 @ sm 6 @ * sm 5 @ sp 6 @ * + 0.5 * + sz 5 @ -0.3 * +
B 7 8 sz 7 @ sz 8 @ * sp 7 @ sm 8 @ * sm 7 @ sp 8 @ * + 0.5 * + sz 7 @ -0.3 * +
B 9 10 sz 9 @ sz 10 @ * sp 9 @ sm 10 @ * sm 9 @ sp 10 @ * + 0.5 * + sz 9 @ -0.3 * +
B 11 12 sz 11 @ sz 12 @ * sp 11 @ sm 12 @ * sm 11 @ sp 12 @ * + 0.5 * + sz 11 @ -0.3 * +
B 13 14 sz 13 @ sz 14 @ * sp 13 @ sm 14 @ * sm 13 @ sp 14 @ * + 0.5 * + sz 13 @ -0.3 * +
B 15 16 sz 15 @ sz 16 @ * sp 15 @ sm 16 @ * sm 15 @ sp 16 @ * + 0.5 * + sz 15 @ -0.3 * +
B 17 18 sz 17 @ sz 18 @ * sp 17 @ sm 18 @ * sm 17 @ sp 18 @ * + 0.5 * + sz 17 @ -0.3 * +
B 19 20 sz 19 @ sz 20 @ * sp 19 @ sm 20 @ * sm 19 @ sp 20 @ * + 0.5 * + sz 19 @ -0.3 * +
B 21 22 sz 21 @ sz 22 @ * sp 21 @ sm 22 @ * sm 21 @ sp 22 @ * + 0.5 * + sz 21 @ -0.3 * +
B 23 24 sz 23 @ sz 24 @ * sp 23 @ sm 24 @ * sm 23 @ sp 24 @ * + 0.5 * + sz 23 @ -0.3 * +
B 25 26 sz 25 @ sz 26 @ * sp 25 @ sm 26 @ * sm 25 @ sp 26 @ * + 0.5 * + sz 25 @ -0.3 * +
B 27 28 sz 27 @ sz 28 @ * sp 27 @ sm 28 @ * sm 27 @ sp 28 @ * + 0.5 * + sz 27 @ -0.3 * +
B 29 30 sz 29 @ sz 30 @ * sp 29 @ sm 30 @ * sm 29 @ sp 30 @ * + 0.5 * + sz 29 @ -0.3 * +
B 31 32 sz 31 @ sz 32 @ * sp 31 @ sm 32 @ * sm 31 @ sp 32 @ * + 0.5 * + sz 31 @ -0.3 * +
B 33 34 sz 33 @ sz 34 @ * sp 33 @ sm 34 @ * sm 33 @ sp 34 @ * + 0.5 * + sz 33 @ -0.3 * +
B 35 36 sz 35 @ sz 36 @ * sp 35 @ sm 36 @ * sm 35 @ sp 36 @ * + 0.5 * + sz 35 @ -0.3 * +
B 37 38 sz 37 @ sz 38 @ * sp 37 @ sm 38 @ * sm 37 @ sp 38 @ * + 0.5 * + sz 37 @ -0.3 * +
B 39 40 sz 39 @ sz 40 @ * sp 39 @ sm 40 @ * sm 39 @ sp 40 @ * + 0.5 * + sz 39 @ -0.3 * +
B 41 42 sz 41 @ sz 42 @ * sp 41 @ sm 42 @ * sm 41 @ sp 42 @ * + 0.5 * + sz 41 @ -0.3 * +
B 43 44 sz 43 @ sz 44 @ * sp 43 @ sm 44 @ * sm 43 @ sp 44 @ * + 0.5 * + sz 43 @ -0.3 * +
B 45 46 sz 45 @ sz 46 @ * sp 45 @ sm 46 @ * sm 45 @ sp 46 @ * + 0.5 * + sz 45 @ -0.3 * +
B 47 48 sz 47 @ sz 48 @ * sp 47 @ sm 48 @ * sm 47 @ sp 48 @ * + 0.5 * + sz 47 @ -0.3 * +
B 49 50 sz 49 @ sz 50 @ * sp 49 @ sm 50 @ * sm 49 @ sp 50 @ * + 0.5 * + sz 49 @ -0.3 * +
B 51 52 sz 51 @ sz 52 @ * sp 51 @ sm 52 @ * sm 51 @ sp 52 @ * + 0.5 * + sz 51 @ -0.3 * +
B 53 54 sz 53 @ sz 54 @ * sp 53 @ sm 54 @ * sm 53 @ sp 54 @ * + 0.5 * + sz 53 @ -0.3 * +
B 55 56 sz 55 @ sz 56 @ * sp 55 @ sm 56 @ * sm 55 @ sp 56 @ * + 0.5 * + sz 55 @ -0.3 * +
B 57 58 sz 57 @ sz 58 @ * sp 57 @ sm 58 @ * sm 57 @ sp 58 @ * + 0.5 * + sz 57 @ -0.3 * +
B 59 60 sz 59 @ sz 60 @ * sp 59 @ sm 60 @ * sm 59 @ sp 60 @ * + 0.5 * + sz 59 @ -0.3 * +
B 61 62 sz 61 @ sz 62 @ * sp 61 @ sm 62 @ * sm 61 @ sp 62 @ * + 0.5 * + sz 61 @ -0.3 * +
B 63 64 sz 63 @ sz 64 @ * sp 63 @ sm 64 @ * sm 63 @ sp 64 @ * + 0.5 * + sz 63 @ -0.3 * +
B 63 64 sz 64 @ 0.3 *

Next, apply the initial excitation to our ground state and truncate the created state to a discarded weight of 10-12 (corresponding to an error of approximately \(L \cdot \sqrt{10^{-12}} = 6.4 \cdot 10^{-5}\):

737: $ syten-apply-op -i gs -o exc -l sp:'sm 32 @'
738: $ syten-truncate -w 1e-12 -i exc -o exc-trunc

and start the time evolution:

739: $ syten-tebd-swappable -l sp -i exc-trunc -s tebd-swappable-file -o tebd -d 0.05 -T 20 --symmetrise true \
740:   -w 1e-12 --apply-method itrunc

This will create a series of .state files. To evaluate the effect of our time evolution, we will use a little Bash snippet to evaluate the local expectation values for each step. Note that to get the correlator, we will also have to multiply the value at each step by the phase factor coming from the second time evolution, viz. \(e^{i \hat H t}\), which is not done here:

741: for i in $(seq -f %g 0.05 0.05 20); do
742:   while [ ! -e tebd_T-${i}.state ]; do
743:     sleep 1;
744:   done;
745:   syten-localexpectation -qs gs sp sp tebd_T-${i}.state > tmp;
746:   printf -- "$(awk '{print $2}' tmp | tr "\n" " ")\n" >> real;
747:   printf -- "$(awk '{print $2}' tmp | tr "\n" " ")\n" >> imag;
748: done

The resulting files real and imag containing the real and imaginary parts of the correlator, one timeslice per line and one site per column, can be plotted using gnuplot, for example:

749: $ gnuplot
750: 
751:         G N U P L O T
752:         Version 5.0 patchlevel 3    last modified 2016-02-21 
753: 
754:         Copyright (C) 1986-1993, 1998, 2004, 2007-2016
755:         Thomas Williams, Colin Kelley and many others
756: 
757:         gnuplot home:     http://www.gnuplot.info
758:         faq, bugs, etc:   type "help FAQ"
759:         immediate help:   type "help"  (plot window: hit 'h')
760: 
761: Terminal type set to 'qt'
762: gnuplot> set terminal pdfcairo  transparent enhanced font "Georgia, 10" fontscale 0.5 size 28cm, 20cm  linewidth 2
763: Terminal type set to 'pdfcairo'
764: Options are ' transparent enhanced font "Georgia, 10" fontscale 0.5 size 28.00cm, 20.00cm  linewidth 2'
765: gnuplot> set output "tebd-local.pdf"
766: gnuplot> set xlabel "Position"
767: gnuplot> set ylabel "Time"
768: gnuplot> set yrange [0:20]
769: gnuplot> set xrange [-32:32]
770: gnuplot> plot "real" matrix u ($1-32):($2/20):3 with image
771: gnuplot> plot "imag" matrix u ($1-32):($2/20):3 with image
772: gnuplot> set logscale cb
773: gnuplot> set cbrange [1e-3:1]
774: gnuplot> plot "real" matrix u ($1-32):($2/20):(abs($3)) with image
775: gnuplot> plot "imag" matrix u ($1-32):($2/20):(abs($3)) with image
776: gnuplot> set output
777: gnuplot> quit

The real part and imaginary part actually look very similar, so we will only show each once:

real-spsm.png

imag-spsm-log.png

To actually calculate the dynamical structure factor, further steps would be necessary:

  • one would have to handle the phase factor arising from the second evolution operator \(e^{i \hat H t}\) to the left of \(\hat S^{+}_i\)
  • not only \(\langle 0 |\hat S^+_{i}(t) \hat S^{-}_{32}(0)|0\rangle\) but also \(\langle 0|\hat S^{-}_{i}(t) \hat S^{+}_{32}(0)|0\rangle\) needs to be evaluated and the results added together.
  • Fourier transformations into momentum space (usually unproblematic, if the system is large enough) and frequency space (problematic, since the correlators usually do not decay within the achievable time frame) need to be implemented.

5.1 Entanglement growth

During the time evolution, you may have noticed a slow-down later on. Additional, the state files have grown even though we have applied the same truncation criterion at all steps:

778: $ ls -1sh tebd_T-?.state tebd_T-??.state | sort -h -k 1
779:  188K tebd_T-1.state
780:  192K tebd_T-2.state
781:  216K tebd_T-3.state
782:  264K tebd_T-4.state
783:  344K tebd_T-5.state
784:  448K tebd_T-6.state
785:  600K tebd_T-7.state
786:  776K tebd_T-8.state
787: 1004K tebd_T-9.state
788:  1.3M tebd_T-10.state
789:  1.6M tebd_T-11.state
790:  2.2M tebd_T-12.state
791:  2.7M tebd_T-13.state
792:  3.2M tebd_T-14.state
793:  3.8M tebd_T-15.state
794:  4.6M tebd_T-16.state
795:  5.4M tebd_T-17.state
796:  6.4M tebd_T-18.state
797:  7.4M tebd_T-19.state
798:  8.5M tebd_T-20.state

This is a symptom of the aforementioned, problematic, entanglement growth. We can quantify it further by calculating the Shannon entanglement entropy at the central bond at each time step:

799: $ for i in $(seq -f %g 0.05 0.05 20); do printf "${i} $(syten-info -qe 1 tebd_T-${i}.state | grep '^  *32 ' | awk '{print $2}')\n"; done

which will show an initial minor drop but the quickly double from initially approx. 0.7 to 1.3. With the spreading excitation, this entropy will grow at every bond of the system, as can be visualised again:

800: for i in $(seq -f %g 0.05 0.05 20); do
801:     while [ ! -e tebd_T-${i}.state ]; do
802:         sleep 1;
803:     done;
804:     syten-info -qe 1 tebd_T-${i}.state | awk '{print $2}' | tr "\n" " " >> entropy;
805:     printf "\n" >> entropy;
806: done

(of course, once the excitation reaches the edge of the system, we would require a larger system if we wanted to evolve physically meaningfully further.) Graphically, the entanglement cone should look like this:

entropy-scaling-local-quench.png

While here our entanglement growth started in the centre and spread through the system, there are cases where it will grow everywhere equally (and hence much worse for the overall simulation): global quenches.

6 Time evolution: global quench and entanglement growth

Briefly, a global quench is a situation where the system reacts in its entirety to a new Hamiltonian governing its dynamics. Examples could be preparation of the ground state of one Hamiltonian (e.g. without a magnetic field) and then a sudden switch-on of that field or simply applying any Hamiltonian as the time-evolution of a random initial state. In the following, we will prepare a new ground state with a smaller magnetic field \(h = 0.05\) and then evolve the system with the same Hamiltonian as in the previous section, still specified in the tebd-swappable-file file:

807: $ syten-dmrg -l sp:'Hp Hz + Hz_s 0.05 * +' -i rnd -s "(m 200 t 0)" -o dmrg-0.05 -f gs-0.05
808: $ syten-truncate -w 1e-12 -i gs-0.05 -o gs-0.05-trunc
809: $ syten-tebd-swappable -l sp -i gs-0.05-trunc -s tebd-swappable-file -o tebd-global -d 0.05 -T 1 --symmetrise true \
810:   -w 1e-12 --apply-method itrunc

We can evaluate the Shannon entropy at every bond exactly the same as before:

811: for i in $(seq -f %g 0.05 0.05 20); do
812:     while [ ! -e tebd_T-${i}.state ]; do
813:         sleep 1;
814:     done;
815:     syten-info -qe 1 tebd-global_T-${i}.state | awk '{print $2}' | tr "\n" " " >> entropy-global;
816:     printf "\n" >> entropy-global;
817: done

You will notice that the time evolution slows down substantially when reaching times \(t \approx 6\ldots 7\). Inspecting just the file sizes of the generated states confirms our suspicion:

818: $ ls -1sh *T-?.state | sort -h
819: 692K tebd_T-1.state
820: 860K tebd_T-2.state
821: 1.6M tebd_T-3.state
822: 3.3M tebd_T-4.state
823: 6.1M tebd_T-5.state
824:  12M tebd_T-6.state
825:  22M tebd_T-7.state
826:  40M tebd_T-8.state
827:  70M tebd_T-9.state
828: 121M tebd_T-10.state

For this reason, I've also limited the generated entropy plot to the range \(t \in [0, 10]\):

entropy-global.png

As can be seen, the entropy grows nearly homogeneously throughout the system, though locally slightly slower than in the case of the local quench.

7 Bosonic systems and local basis optimisation [extra]

For many one-dimensional systems, the linear geometry of matrix-product states are quite perfectly suitable. One class of cases where a different choice of geometry makes a major difference are bosonic systems with many local degrees of freedom (e.g. up to \(n = 50\) bosons per site). First, create the lattice and an initial random state (the syten-mps-boson-spinless-nil lattice generator is not very optimised and we are not using any symmetries now, so this might take a while):

829: $ SYTEN_DEFAULT_TPO_TRUNCATE=p syten-mps-boson-spinless-nil -l 20 -n 50 -o lattice
830: $ syten-random -l lattice -o rnd.mps

We will now calculate the ground state of this system with local on-site interaction strength 0.1:

831: $ SYTEN_DEFAULT_TPO_TRUNCATE=p syten-dmrg -l lattice:'Hu 0.1 * Hj +' -i rnd.mps \
832:   -s "(m 10 t 0 x 5) (m 20) (m 40) (m 60) (m 80) (m 100) (m 120) (m 140) (m 200)" -o dmrg.mps

The calculation will take a few minutes, so open a new window and generate a local-basis-optimised state of the initial random state (see here or or here) and start another DMRG run, this time also specifying a ‘secondary’ bond dimension with sm and a ‘secondary’ truncation threshold with ‘st’:

833: $ syten-mps2lbo -i rnd.mps -o rnd.lbo
834: $ SYTEN_DEFAULT_TPO_TRUNCATE=p syten-dmrg -l lattice:'Hu 0.1 * Hj +' -i rnd.lbo \
835:   -s "(m 10 t 0 x 5 sm 5 st 0) (m 20 sm 10) (m 40 sm 20) (m 60) (m 80) (m 100) (m 120) (m 140) (m 200)" -o dmrg-lbo

You will notice that the LBO-MPS-DMRG calculation relatively quickly overtakes the already-running MPS-DMRG calculation and will print lower energies (in the eigth and eleventh column). The size of its generated states will also be considerably smaller (some 50%).

Since the entire system most likely needs a bond dimension larger than 150, which, without symmetries, is considerable numerical effort, you should probably stop both calculations now using Ctrl-C.

Author: Claudius Hubig

Email: c.hubig@physik.uni-muenchen.de

Created: 2017-08-24 Do 15:00

Emacs 24.5.1 (Org mode 8.2.10)

Validate