{'fraction_non_alphanumeric': 0.034078127374981, 'fraction_numerical': 0.008511931904544764, 'mean_word_length': 4.603134048714018, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Matrix Product State applications for the ALPS project

The density-matrix renormalization group method has become a standard computational approach to the low-energy physics as well as dynamics of low-dimensional quantum systems. In this paper, we present a new set of applications, available as part of the ALPS package, that provide an efficient and flexible implementation of these methods based on a matrix-product state (MPS) representation. Our applications implement, within the same framework, algorithms to variationally find the ground state and low-lying excited states as well as simulate the time evolution of arbitrary one-dimensional and two-dimensional models. Implementing the conservation of quantum numbers for generic Abelian symmetries, we achieve performance competitive with the best codes in the community. Example results are provided for (i) a model of itinerant fermions in one dimension and (ii) a model of quantum magnetism.


Introduction
The density matrix renormalization group method (DMRG) [3] is currently the method of choice to investigate low-dimensional systems in condensed matter physics. Its applications range from the study of superconductivity, quantum magnetism and exotic phases of matter to simulation of quantum circuits.
On the theoretical side, a lot of progress has been made in understanding the success of this method on the basis of matrix product states (MPS), which is the class of variational states underlying the DMRG method. This has facilitated the development of new algorithms, e.g. for real-time evolution [4][5][6], mixed-state (finite temperature) calculations [7,8] as well as their combination [9][10][11]. These have been formulated directly in the language of MPS, and are now widely being applied and further developed.
Recently, there has also been a surge of applications of DMRG to two-dimensional systems [12]. While it was realized early on that DMRG is exponentially hard for systems in more than one dimension [13], modern DMRG codes have become sufficiently efficient to combat this exponential scaling and reach system sizes where reliable conclusions about the thermodynamic limit of the 2d system become feasible.
In this paper we report the release into the ALPS project [1,2] of a set of MPS codes developed within the Swiss HP2C initiative [14], whose target was to develop new efficient and productive applications for current and future high performance computers. These applications have been designed to achieve the high performance required for the study of two-dimensional systems, and additionally to bind to the generic and well-known ALPS model and lattice libraries, which makes our MPS codes appealing for researchers in various fields. The new code has thus far been employed to study condensed matter physics problems such as ordering in an SU(3) Heisenberg model [15], the analysis of critical points for a system of supersymmetric fermions [16], the analysis of interacting mesoscopic structures [17], spin liquid phases in a Kagomelattice spin model [18], and recent method developments such as the applications of DMRG as impurity solver within DMFT [19] and the development a multigrid DMRG scheme [20], A number of open-source implementations of DMRG are currently available [1,[21][22][23][24]. With this release, we add a set of applications that, while implementing a generic framework amenable to many extensions, offers performance on par with the best available codes. Furthermore, it is seamlessly integrated with the widely used ALPS libraries.
The outline of this paper is as follows: In Section 2 we first focus on the characteristics and conventions of the new MPS code. In Section 3 we explain the main application parameters and show the results of a tutorial. For the theoretical background on MPS and DMRG we point the reader to the review [25] and references therein.

Algorithms and implementation
The matrix product state (MPS) is a variational ansatz wave function which for a chain of N sites reads |ψ = σ where the exponentially growing number of parameters of the full quantum state is efficiently truncated to a polynomial number ∝ LM 2 contained in the M × M matrices A σi i . This introduces an approximation that can be controlled by the matrix size M . Matrix-product states are an example of a tensor network state. Such states are compactly represented graphically, as shown in Fig. 1. A tensor network is represented as a graph where vertices represent tensors, and edges represent indices that are contracted. Thus, in the upper panel of Fig. 1, the blue circles represent the tensors (A σi i ) αi−1αi of Eqn. (1); the edges connecting the tensors represent the indices that are summed over to yield matrix multiplication; and the open vertices pointing down represent the physical indices σ i . The lower panel of the figure shows a closed tensor network, i.e. without any open legs, that thus corresponds to a scalar. The contraction shown yields the calculation of an expectation value for an operator represented by a matrix-product operator (see Section 2.1).
It has been shown that ground states of gapped Hamiltonians in one dimension are accurately represented by a matrix-product state with a matrix size that grows only polynomially in the system size [26,27]; if only local properties are required, one can obtain even better scaling [28]. Intuitively, this result rests on the fact that an area law holds for these systems [29], i.e. the bipartite entanglement entropy in these states is a function only of the area of the cut between two regions and not the total system size. In one dimension, this leads to constant bipartite entanglement entropy for gapped states.
The MPS ansatz is inherently a one-dimensional open chain. To simulate systems on periodic systems, ladders or even two-dimensional systems, a (a)Î 1Î2ĉ † 3ĉ 4 mapping from these lattices to a chain has to be chosen. This will lead to a Hamiltonian with longerranged interactions on the chain, where the range of interactions is related to the width of the ladder or 2d system. More importantly, the entanglement structure of the state is affected. The area of a cut bipartitioning the system into two halves will generally grow with the linear size L of the system. Since the bond dimension required to capture entropy S in an MPS is exponential, M ∼ e S , this implies an exponential cost for MPS simulations for twodimensional systems. The specific mapping from the 2d system to the chain therefore may have great impact on the accuracy and performance of simulations, and should therefore be chosen very carefully to minimize entanglement across each cut.

Representation of Matrix Product Operators
To introduce the notion of matrix product operators, it is useful to consider a 1d quantum lattice system with the same Hilbert space H on each site. Let {Ô β k } be the set of Hermitian operators that can act on the k'th site; as an example, these are identity matrix and the three Pauli matrices, {I, σ x , σ y , σ z } for the case of S = 1/2 spins. We can then write any Hermitian operatorQ acting on the entire lattice aŝ This makes use of the fact that the operators themselves form a Hilbert space, analogous to quantum state. The above expression suggests that the coefficients c(α 1 , . . .) can be decomposed into a matrix product, analogous to Eqn. 1. This representation is referred to as matrix-product operator (MPO). Analogous to a matrix-product state, such a matrix-product operator has an associated matrix size. Local Hamiltonians can generally be represented as an MPO of low bond dimension; even for long-range Hamiltonians, an efficient MPO description is often possible. Consider now writing a local Hamiltonian, say for free fermions, in this form. For simplicity, let's consider this model on three sites but with additional nextnearest neighbor hopping: The most naive way to perform this decomposition is to assign an independent auxiliary index to each term in each of the sums: This leads to an MPO of bond dimension 6, which is very inefficient and redundant. A more efficient version is the following: This MPO implements the same Hamiltonian, but has only bond dimension 3 by making use of the fact that the operatorsĉ † 1 ,ĉ 1 as well asĉ † 3 ,ĉ 3 occur both in the nearest-as well as next-nearest neighbor hopping term. Note that in the notation above, we have taken the elements of the matrices to be fermionic operators themselves. The antisymmetry of the central matrix in this case is a special property of the 3-site ring and reflects the anticommutation of fermionic operators. In the actual implementation, fermionic statistics is implemented via a Jordan-Wigner transformation.
Our code implements an optimized MPO network that merges common branches, as illustrated in the above example, thus reducing the matrix size of the MPO substantially. This is sketched in the lower panel of Figure 2 (for a more extensive explanation of similar notation, see Ref. [30]). Note that in principle, the same compression techniques that are available for reducing the matrix size of MPSs can also be applied to MPOs; however, in some cases constructing the MPO in a naive way and then compressing is prohibitively time-consuming. Therefore, the code is written to perform many optimizations when the MPO is first constructed.

Fermionic operators
Special treatment is needed to account for the anti-commutation relations {ĉ(i),ĉ † (j)} = δ ij characteristic of fermions. All fermionic operators have to be converted into a string of bosonic operators according to the Jordan-Wigner transformation [31]:ĉ wheref (k) is the fermion sign operator, a diagonal matrix with a −1 entry for each local basis state with an odd number of fermions, and 1 for each local basis state with an even (including 0) number of fermions, andb(i) are hard-core boson operators satisfying the commutation relation [b(i),b † (j)] = δ ij and (b(i)) 2 = 0. The fermion sign matrixf (k) can be written asf (k) = 1 − 2b † (k)b(k). To represent spinful fermions, two flavors of bosonsb 1 (k),b 2 (k) have to be introduced and a normal-ordering of operators has to be chosen; in the ALPS libraries, the ordering is chosen aŝ b 1 (1)b 2 (1)b 1 (2)b 2 (2).... The Jordan-Wigner string must then include fermion sign operators for both flavors of bosons, i.e. two per physical site. Noting thatf 2 = I and [f (i),b ( †) (j)] = 0 for i = j, we find that for a fermionic hopping term the transformation simplifies to just a sequence of f (k) matrices between the two sites i and j: A similar simplification will apply to other types of fermionic bond terms; bond terms where each term is diagonal, such as density-density interactions, do not require Jordan-Wigner transformation.

Variational optimization of the MPS
Given the variational nature of the MPS ansatz, the ground state of a model described by the Hamil-tonianĤ is obtained by an optimization process of all variational parameters looking for the energy minimum.
The most common way to optimize an MPS, as already described by S. White in '92 [3] is to optimize one or two tensors at the time, keeping all the others constant. For the case of the single site optimization the equation to solve is which maps to solving for the lowest eigenvector of the generalized eigenvalue problem MA i −λN A i = 0, which is solved iteratively with the Jacobi-Davison algorithm provided by the IETL package in ALPS. Note that in the iterative solver the matrix M is never computed, but we efficiently contract all tensor legs on the fly in order to keep a complexity O(M 3 ). Convergence is reached after sweeping a few times through the system. One sweep is here defined as optimizing twice all tensors in the chain, i.e. moving from left to right and then from right to left.
In our mps_optim application we implement the so called single-site and two-site optimizations. Figure 3: Construction of the environment tensor used as a local orthogonality constraint. The orthogonality constraint shown in the left network corresponds to a mixed MPS tensor which is used to evaluate ψ 0 |ψ 1 in the local optimization.
The first one optimizes one tensor at the time, whereas the latter one fuses two sites together to optimize the bond between them. The user is free to choose between them by changing the optimization input parameter to singlesite or twosite, respectively. It is strongly advised to always play around with both, in order to find what fits better the system to be solved.
In the case of the single site optimization we implement the procedure proposed in Ref. [32]. Here, the two shortcomings of the single-site optimization -the possibility of trapping on local minima and the inability to change the bond dimension dynamically -are remedied by enlarging the reduced density matrix after each optimization using Assuming a partition of the system into a sub-block A and sub-block B such that the Hamiltonian is b l , the first term inρ A is the usual reduced density matrix of the subsystem A, whereas the second term introduces some reshuffling of the truncated states that protects against trapping in the optimization into local minima. The parameter α is a small number -typically 10 −4 -which is slowly taken to zero during sweeping. In our implementation we allow three values for it alpha initial, alpha main and alpha final (see the parameters in Table 2 for further details).

Excited states
In traditional DMRG approaches, the calculation of excited states is implemented through multi-state targetting, where the reduced density matrix after each iteration is truncated to describe a number of low-lying states simultaneously. While this approach has been shown to work reasonably, a more controlled approach is to sequentially find low-lying states by first finding the ground state |ψ 0 , then finding the lowest-energy state |ψ 1 that is orthog-onal to the ground state, ψ 0 |ψ 1 = 0, etc. This approach was first described in Ref. [33].
In the MPS context, the overlap of two states can be efficiently calculated; furthermore, when optimizing a local tensor, the orthogonality constraint can be expressed as a local constraint of orthogonalizing the MPS tensor against the environment tensor shown in Fig. 3. For higher excited states, several such constraints need to be taken into account. In an iterative eigensolver, such as the Jacobi-Davidson method used here, these constraints can be taken into account by projecting newly generated Krylov vectors (or, in the case of Jacobi-Davidson, search vectors) into the orthogonal complement of the space spanned by the constraints. When sweeping through the state, the environment of Fig. 3 can be updated at very low computational cost.

Time evolution
To simulate time evolution, we follow the approach of Refs. [5,6,34]. The key insight to performing time evolution using MPS is that instead of time evolution, one can also variationally optimize the difference between an MPS |ψ 1 and an MPS with a unitary operator applied, U |ψ 0 ; that is, the optimization problem can be solved efficiently over normalized MPS.
Here, the operator U can either act locally on a few sites, or be a sum of terms expressed as an MPO. Another key principle is that the evolution under a local Hamiltonian,Ĥ = nĥ n , can be expressed as a product over local unitaries using a Trotter-Suzuki decomposition [35,36]. Here, the time evolution operator exp(−iĤt) for a small time step ∆t is decomposed into multiple products of the non-communing terms in the Hamiltonian. To first order the Trotter-Suzuki decomposition is (15) whereĥ n can themselves be sums of terms that do not overlap and hence commute. In the case of a simple nearest-neighbor Hamiltonian on a chain,ĥ 1 andĥ 2 are the terms acting on even and odd bonds, respectively. In general better accuracy is obtained by second and fourth orders which are implemented in the mps_evolve application (see Appendix A.1).
Naively, the application of the unitary to an MPS increases the bond dimension to M = M · D h , where D h > 1 depends on the structure of the Hamiltonian; for details, see below. Therefore, after each contraction, the MPS matrices are truncated to their original size either with an singular value decomposition, or with a variational approximation (see Ref. [25]). Our application implements two types of tensor contractions for time evolution, nearest neighbors and mpo.
Nearest neighbors. Nearest neighbors time evolution algorithm is an optimized version used when the Hamiltonian contains only nearest neighbors bond term (in the case of ladders or 2d systems after having been mapped it to a Hamiltonian on a chain). The bond gates are contracted with the MPS tensors and a truncation via singular value decomposition is applied on the fly, following the TEBD method described in Ref. [34].
MPO. MPO time evolution is a more generic algorithm where each exp(−iĥ n ∆t) is transformed into an MPO, which allows to encode bond terms between non-neighboring sites. After multiplying the MPS wave function with the MPO, this is truncated to its original size M with a variational compression. Special care has to be taken for fermionic long-range operators, this procedure is explained in Appendix A.2.

Efficient tensor storage with abelian symmetries
Symmetries manifest themselves in many physical models, i.e. there are global operatorsĜ that commute with the system Hamiltonian, [Ĥ,Ĝ] = 0. In such a case, due to the fact that these operators can be simultaneously diagonalized,Ĥ becomes block-diagonal in a basis of eigenstates of G. The blocks of the Hamiltonian are labelled by quantum numbers, corresponding to possible eigenvalues ofĜ. It is often computationally favorable to diagonalize the Hamiltonian separately within each quantum number sector.
In the case of DMRG, this is particularly easy when the symmetry has local generators, such that the reduced density matrices for any part of the system can also be block-diagonalized. Since the bond indices of the matrix product state essentially enumerate eigenstates of reduced density matrices, In the pair (|k , s) we describe a charge block with label |k (k particles) of size s. The MPS always starts with the trivial block (|0 , 1) and ends with the total charge sector that the user wants to target, n particles in this example.
these can be associated with quantum number sectors. Each matrix of the MPS then obeys the symmetry locally, i.e. there are local conservation laws, and the matrices also become block-diagonal. For each bond of the tensor network, a gauge choice has to be made fixing the direction of the bond. When considering, for example, particle number conservation, this is a choice of which bonds are considered going into and which coming out of the node of the network [37]. The choice we make in our implementation is shown in Fig. 4. The MPS always starts with a trivial index on the left (e.g. only one sector with zero particles), new states are then added on every site from the local Hilbert space σ i , and eventually the last tensor is fixed to the total quantum number sector chosen by the user.
Given our convention, locally the conservation of quantum numbers reads where α i , σ i are quantum numbers and • represents the group product. The most convenient representation for the MPS tensors is in a matrix-like form where the physical index σ i is fused together with one of the auxiliary indexes α i or α i+1 . We then profit from the fact that the matrix A σiαi; αi+1 is block-diagonal in the quantum number labels, and we store only the dense blocks. Linear algebra operations acting on these dense matrices are dispatched to optimized BLAS and LAPACK libraries available on the user machine 1 .
The MPS codes presented in this paper are implemented to make use of Abelian symmetries. Dealing with conserved quantum numbers involves an intense calculation of indexes labels. We therefore set the symmetry at compile-time to allow for a maximum amount of compiler optimization. By default, the code is compiled for U (1) ⊗n symmetry, where the value n is fixed during compilation by the macro DMRG_NUMSYMM, as well as Z 2 symmetry 2 . Other symmetries may easily be added.

Parallelization
We parallelize our applications for shared memory architectures using the OpenMP compiler directives. The parallel for loop directive is employed, for example, in the contraction of the local MPS tensor with the local MPO terms. In this routine several loops over the size of the auxiliary MPO dimension operate on independent parts. Since the MPO auxiliary dimension grows linearly with the width of a 2d system, we profit from more parallelism for large systems. The OpenMP parallelization is enabled at compile-time with the CMake option ALPS_ENABLE_OPENMP, and can be controlled at run-time via standard environment variables such as OMP_NUM_THREADS.
Many vendors provide support for multi-threaded BLAS and LAPACK libraries. For large matrices sizes our application will profit automatically from this parallelization, too. However, special care must be taken when using both the OpenMP parallelization and a parallel BLAS/LAPACK library simultaneously; this will only lead to improved scaling when using a compiler that supports such nested parallelism.

Codes and examples
Our set of MPS tools is composed of an application for the calculation of ground states and low energy excitations (mps_optim), a time evolution application (mps_evolve) and two utilities (mps_meas and mps_overlap) to perform measurements on the wave functions generated with the previous programs.
The interface of the applications is consistent with other ALPS applications, i.e. the same simulation can easily be run with exact diagonalization, the MPS optimization or even Monte Carlo applications. On the outer level a scheduler reads an XML input file listing all parameter sets to be computed.
These parameter files can be created with the standard ALPS tools: (1) converting text files to XML with the parameter2xml command line tool, (2) using the convenient Python functions available in the pyalps package or (3) generating a workflow from VisTrails [38].
Each simulation generates three output files: simulation.out.h5 contains the simulation parameters, the iteration results and the final measurements in HDF5 format (see Appendix C for the detailed HDF5 schemes); simulation.out.xml contains the simulation parameters and the final results in XML format conforming with the ALPS schemas [1,2]; simulation.out.chkp is the simulation checkpoint file containing the simulation status and the final MPS wave function that can be used to restore the same simulation or as input for other calculations.
Together with the applications, in the ancillary files, we ship full example simulations to illustrate how to create parameter files, run simulations and analyze the data. In sections 3.2 and 3.3 we describe some of these results.

Input parameters
In Table 1, Table 2 and Table 3 we list all parameters and the default values used by the applications. They are split between common MPS parameters (valid in all codes), optimization parameters (available only in mps_optim) and time evolution specific parameters (for mps_evolve). Models, lattice and measurements are parameters of the ALPS model and lattice libraries, we point the reader to the ALPS references [1, 2] for further details.

Initial states
The MPS applications can be initialized with four different types of initial states.
default. The default initial state constructs an MPS which contains all allowed quantum number sectors that fulfill the total quantum numbers chosen by the user. In each MPS tensor all sectors will have a maximum size of 5 × 5 and they will be filled with random numbers.
thin. The thin initial state is an extension of the default option. After constructing the default MPS the routine performs a compression retaining only a total 20 truncated states per site. This initial state reduces considerably the number of sectors that the code has to tackle in the warm up sweeps. Sometimes it leads to faster convergence. basis state. A many body basis state is easily encoded in a product state as an MPS with bond dimension M = 1. With this option the user can design its own initial state listing the desired local quantum numbers in the input file as illustrated in Table 1.
input file. The checkpoint file generated by a simulation can easily be reused in a new simulation by specifying it in the initfile parameter. This option is often used to start a time evolution from the ground state of an other model.

Temporary memory
An MPS simulation easily requires more than 100 GB of memory for the storage of boundary objects, i.e. partial contractions of the MPS network needed to reduce the complexity of the variational optimization to only O(M 3 ).
Our codes are capable of efficiently move this temporary memory to disk. This option is enabled by giving a non-empty value to the parameter storagedir. Since the application will often write > 100GB of temporary data it is advised to use fast filesystems, e.g. running on a local hard disk instead of a network mount. Once a simulation finishes the user can safely erase all temporary data.
whereĉ † ↑ (i) andĉ † ↓ (i) create a fermion at site i with spin up and down, respectively, andn σ (i) = c † σ (i)ĉ σ (i) counts the number of particles with spin σ ∈ {↑, ↓} on site i. For a relatively strong repulsive interaction U/t = 8 and an average filling n = 0.875 this model is expected to fall into the Luther-Emery universality class [39,40] with a spin gap and a gapless charge excitation. To leading order both the long-range pairing correlation function D(|i − j|) = D(l) ∼ l −1/Kρ -defined as with 2)] creating a singlet on the ith-rung -and the long-range charge density wave correlation function C(|i−j|) = C(l) ∼ l −Kρ -defined as -decay algebraically with an exponent that depends only on the non-universal parameter K ρ , which we have to determine from the simulation at the chosen model parameters.
Long-range correlations have to propagate through the whole system and are therefore easily suppressed when truncating to only a few DMRG states. To treat the DMRG truncation correctly we run several instances of the simulation with increasing bond dimension M , until the observables converges or we can confidently extrapolate their value to M → ∞.
In this example we simulate a 2 × 96 ladder with bond dimension varying from a computationally cheap M = 800 (convergence takes ca. 2.3h on 5 cores) to a computationally expensive M = 3600 (convergence takes 46h on 5 cores). The run script is:  Because of the iterative optimization in DMRG, one should always check for convergence with the number of sweeps through the system. In Fig. 5, we see that both the ground state energy and the truncated weight have converged in our simulations. Note that in the above parameters we enable the calculation of the energy variance Var[Ĥ] = Ĥ 2 − Ĥ 2 . This is a very convenient quantity for extrapolating observables to M → ∞, corresponding to Var[Ĥ] → 0 for an eigenstate ofĤ. As an example the ground state energy depends to first order linearly on Var[Ĥ]. In Figure 6 we show the extrapolation of the energy generated with the following code: The evaluation of the pairing correlation function 2) which was listed as a correlation measurement in the parameters listed above. Since field_du and fielddag_ud are defined in the file mymodel.xml 3 as bond operators, the application applies the four operators to all pairs of nearest neighbors sites, i.e. i, i + 1, j, j + i for all i = 1, . . . , N − 4 and j = i + 2, . . . , N , effectively scaling as O(N 2 ). Form the measured D(i, j) we compute D(l) = D(|i − j|) by averaging over a few rungs at distance l = |i−j| around the center of the ladder. This technique averages out oscillations induced by the open boundaries. Results are shown in Figure 7, where one notes a very slow convergence of the pairing correlation function D(l) with the retained DMRG states; for M = 800 pairing at large distances is underestimated. Large enough MPS matrices, like the M = 3600 data, are needed to correctly analyze the system.
From this example simulation we observe that D(l) decays with an exponent consistent with K ρ > 1, which is the condition for a dominant superconducting phase in the Luther-Emery universality class.

Example 2: Evolution of a domain wall
As a second example we run a simple time evolution simulation of the Heisenberg spin XX model, described by the Hamiltonian with S x (i), S y (i) the x, y component of the spin−1/2 operator S(i).
The time evolution of an initial domain wall |ψ(t = 0) = |↓↓ . . . ↓↑ . . . ↑↑ -a state with all up spins on the right side of the chain and the down spins on the left -under the spin XX Hamiltonian was solved exactly by a Jordan-Wigner transformation to free fermions in Ref. [41]. In particular, the local magnetization at distance n from the center after time t is described by a sum of Bessel functions as with j i (t) is the Bessel function of order i. The setup script is very similar to the one previously used for the ground state calculations, for a simulation of 500 time steps up to a total time τ = 20 /J it reads: Note the definition of the initial state initial local Sz listing the value of the S z quantum number on each sites. In Figure 8 the local magnetization calculated with the above simulation is compared with the exact solution from Eqn. 20 showing very good agreement. To produce a similar figure one simply has to load the results for each time step with the loadIterationMeasurements function and select the magnetization at distance n from the center, as shown in the following snippet. A slightly more complex version generating a movie of the local magnetization per site is also provided in the source code of the example.

Acknowledgments
We gratefully acknowledge support by the wider ALPS community. We are indebted to Philippe Corboz and Ulrich Schollwöck for useful discussions in the design process, and to U. Schollwöck for providing a draft of his Python code. We also thank Jan Gukelberger, Hiroshi Shinaoka and Lei Wang for being early users of the framework and applications. Simulations were performed on the PASC Mönch cluster at ETH Zurich. Listing 1: Superconductor with Z 2 symmetry.

/parameters/PARAMETERNAME
Value of the parameter PARAMETERNAME .

Measurements
/spectrum/results/OBSNAME/mean/value Vector of expectation values of the observable OBSNAME , each entry in the vector is a different eigenstate.
/spectrum/results/OBSNAME/labels Labels for the x-coordinate of the observable.

Iterations
/spectrum/iteration/NUM/parameters/ PARAMETERNAME Value of the parameter PARAMETERNAME specific to the iteration NUM . Only modified parameters appear in this group.

/spectrum/iteration/NUM/results/OBSNAME/ mean/value
Expectation values of the observable OBSNAME at the iteration NUM .
/spectrum/iteration/NUM/results/OBSNAME/ labels Labels for the x-coordinate of the observable.
Helper function: pyalps.loadIterations(files) Defines a new average measurement called NAME. Its value specifies the operator to be measured. Example:

MEASURE AVERAGE[Density] = "n" MEASURE CORRELATIONS[NAME ]
Defines a new correlation measurement called NAME. Its value specifies the operators to be correlated separated by ":" (colon) If the operators define two or more BONDOPERATORS, e.g., Op_A and Op_B the correlation measurement compute all observables ψ| Op A (i, i+1) Op B (j, j + 1) |ψ . Example: MEASURE CORRELATIONS[Onebody Density Matrix] = "bdag:b" MEASURE HALF CORRELATIONS[NAME ] Same as MEASURE CORRELATIONS but it does not exchange the order of operators. If the input is, e.g., bdag:b, the first operator will be evaluated at all locatios i ∈ [0, L − 2] but the second operator only at locations j ∈ [i, L − 1].

MEASURE[EnergyVariance]
False Path to the directory to be used for temporary storage. If empty, temporary storage is disabled. Single site only. Number of initial sweeps where the correction factor has the value alpha initial. nmainsweeps Single site only. Number of sweeps after ngrowsweeps where the correction factor has the value alpha main. alpha initial 1e-2 Correction factor for single site optimization [32] in the initial part alpha main 1e-4 Correction factor for single site optimization [32] in the main part alpha final 1e-8 Correction factor for single site optimization [32] in the final part ietl jcd toll 1e-8 Convergence tolerance of the Jacobi-Davidson solver ietl jcd gmres 0 Convergence tolerance of the Jacobi-Davidson solver ietl jcd maxiter 8 Maximum number of iterations in the Jacobi-Davidson solver Values assigned to the parameter P at each time step