# Simulation Details¶

## Including the required Python libraries¶

Python libraries are included in a file by placing

```
import PythonLibrary as pl
```

at the beginning of the file. Here, `PythonLibrary`

is the name of the
Python library to be included, and `pl`

is the
namespace
where the contents of this module are placed. As an example, if we wanted to
include MPSPyLib and call the function `Ops.BuildSpinOperators()`

from
this library, we would use the code

```
import MPSPyLib as mps
Operators = mps.BuildSpinOperators()
```

Note the use of the namespace identifier `mps`

to specify that
`Ops.BuildSpinOperators()`

comes from `MPSPyLib`

. In addition to
MPSPyLib, which is always required for MPS simulations, we will also
frequently include NumPy, as it provides the matrix functionality used
internally for the code. We will use the namespace identifier `np`

for NumPy.

## Specifying the parameters of a simulation¶

A simulation is specified in MPSPyLib as a Python dictionary of simulation data.
A Python dictionary, also called a
hash table, is
an unordered collection of `key:value`

pairs. For example, we could
construct a dictionary which has a single key `'name'`

which maps to the
string `'fred'`

as

```
mydict = {'name' : 'fred'}
```

A dictionary is specified by curly braces, , with
*key:value* pairs separated by a colon, `:`

, which assigns each key to a
value. In python dictionaries, the keys can be any type of immutable data
(i.e. strings, tuples, or numbers). Each key must be unique to the dictionary.
The values assigned to each key is not necessarily unique, and can be any data
type. The value corresponding to a key is accessed by calling that element of
the dictionary using square brackets, , referring
to the dictionary name `mydict`

and the key `'name'`

with the command
`mydict['name']`

.

We will refer to the Python dictionary specifying a simulation in MPSPyLib as
`parameters`

. Every simulation must have the non-optional keys listed in
the table in `tools.WriteFiles()`

. Furthermore optional parameters
are available if the default values are not suitable for the goal of the
simulation.
The particular choice of these parameters defines which algorithms are to be
used and their convergence properties. In addition, `parameters`

should also
contain the values of the Hamiltonian parameters used to define the
Hamiltonian, see Sec. Defining the Hamiltonian and other many-body operators. The simulation specification `parameters`

is converted into simulation data which is readable by the Fortran backend using
the function

```
MainFiles = mps.WriteFiles(parameters, Operators, H)
```

where `Operators`

specifies the operators defining the on-site Hilbert space,
see Sec. Defining local operators, and `H`

specifies the Hamiltonian, see
Sec. Defining the Hamiltonian and other many-body operators. The specification of simulations `parameters`

can also
be a list of dictionaries, which allows for easy batching and parallelization
of simulations. A full set of examples are included in the *Examples*
subdirectory of OSMPS.

## Defining local operators¶

### Overview of operators in MPSPyLib¶

The way that ordinary, meaning non-symmetry-adapted, operators are stored
in MPSPyLib are as dictionaries, which are unordered collections of
key:value pairs (Symmetry-adapted operators are discussed in Sec.
Abelian symmetries). This means that the set of operators is stored as a
dictionary, call it `Operators`

, and

```
Operators['tag']
```

returns the matrix representation of the operator `'tag'`

as a
NumPy array.
The default operators built by the convenience functions are listed in a
table inside of the function description. The options are spin operators
from `Ops.BuildSpinOperators()`

, bosonic operators defined in
`Ops.BuildBoseOperators()`

, fermionic operators created by the
function `Ops.BuildFermiOperators()`

, or operators for a Bose-Fermi
system `Ops.BuildFermiBoseOperators()`

building all operators present
in the bosonic and fermionic set. The total number operator
with the key `'ntotal'`

is built in addition for the latter one. For more
information on the optional arguments of these functions, one can use the
help functionality of Python as alternative to this documentation, e.g.

```
help(mps.BuildBoseOperators)
```

within the Python interpreter.

New operators can be defined by arithmetic on existing operators, and the new operator is designated by a user-defined key. As an example, consider that one wants to construct the bosonic on-site interaction operator

and call it `'Interaction'`

. This may be achieved as

```
import MPSPyLib as mps
import numpy as np
Operators = mps.BuildBoseOperators(Nmax=2)
Operators['Interaction'] = 0.5 * (np.dot(Operators['nbtotal'], Operators['nbtotal']) - Operators['nbtotal'])
```

Here, `np.dot`

is matrix-matrix multiplication from the NumPy libraries.

### Operators with internal degrees of freedom¶

MPSPyLib contains two different methods for generating operators for particles
with internal degrees of freedom. The first is to include a value for the spin
of the particle in `Ops.BuildBoseOperators()`

or
`Ops.BuildFermiOperators()`

. In that case each operator from the set
of operators, with the exception of the total number and Fermi
phase operators, obtains a spin projection index e.g. `tag_m`

, where `m`

runs from 0 to , with the spin. In addition, the number
operators in each component `n*_m`

and the spin operators `'*sz`

,
`'*splus`

, and `'*sminus`

are built, where `*=b`

for bosons and `*=f`

for fermions.

The second means of producing particles with internal degrees of freedom in
MPSPyLib is through specifying a nonzero value for the `nFlavors`

optional
argument of `Ops.BuildBoseOperators()`

or
`Ops.BuildFermiOperators()`

. In this case each operator from the set
of operators, with the exception of the total number operator,
obtains a flavor projection index e.g. `tag_m`

, where `m`

runs from 0
to `nFlavors-1`

, and the number operators in each component `n*\_m`

are
built, where `*=b`

for bosons and `*=f`

for fermions. If both `spin`

and `nFlavors`

are nonzero, then the operators are indexed as
`'tag\_flavor\_m'`

, where `flavor`

runs from 0 to `nFlavors-1`

and `m`

runs from 0 to , with the spin. In
`Ops.BuildFermiBoseOperators()`

different flavors and spins can be
specified for the bosonic and fermionic species. For more information on this
syntax consult the docstring for `Ops.BuildFermiBoseOperators()`

.

Finally, we note that the matrix representations of the operators may be
viewed by using the `print`

functionality of python. That is, if we
wanted to see the bosonic destruction operator `b`

, we could call

```
import MPSPyLib as mps
import numpy as np
Operators = mps.BuildBoseOperators(Nmax=2)
print(Operators['b'])
```

Likewise, we can print all the operators as

```
for i in Operators:
print(Operators[i])
```

## Defining the Hamiltonian and other many-body operators¶

The related class for these type of operators is `MPO.MPO`

. The
natural operator structure for MPS algorithms is a matrix product operator
(MPO). Operators can be constructed as MPOs using a set of rules which specify
how operators acting on local Hilbert spaces are combined to create an operator
acting on the entire many-body system. Using a canonical form for MPOs, several
operators constructed from these rules can be put together as a single
MPO [WC12]. OSMPS uses this essential idea of creating
operators from a set of predetermined rules to construct the Hamiltonian and
other operators appearing in an MPS simulation.

An operator acting on the many-body state, i.e. the MPO, is an object of the
`MPO.MPO`

class in MPSPyLib. For an instantiation of the MPO class,
named `H`

in the following, is obtained as

```
H = mps.MPO(Operators)
```

Once we have an MPO object, we add terms to it using the
`Ops.AddMPOTerm()`

method. In particular, the rules are

`site`

rule (`mpoterm.SiteTerm`

). Example:

H.AddMPOTerm('site', 'sz', hparam='h_z', weight=1.0)generates , where

`Operators['sz']}`

and`Operators`

was already passed when creating the instance of the MPO. The second argument contains the relevant keys of the operators to be used. Because only the names of the operators are used, this routine handles ordinary and symmetry-adapted operators on the same footing.

`bond`

rule (`mpoterm.BondTerm`

). Examples:

H.AddMPOTerm('bond', ['splus','sminus'], hparam='J_xy', weight=0.5)generates , where denotes all nearest neighbor pairs and ,

H.AddMPOTerm('bond', ['sz', 'sz'], hparam='J_z', weight=1.0)generates , and

H.AddMPOTerm('bond', ['fdagger', 'f'], hparam='t', weight=-1.0, Phase=True)generates , where is a fermionic destruction operator. That is, terms with Fermi phases are cared for with the optional

`Phase`

argument and Hermiticity is automatically enforced by this routine. The Fermi phase is added through a Jordan-Wigner string of the operators in`Operator`

with key`FermiPhase`

. The functions`Ops.BuildFermiOperators()`

and`Ops.BuildFermiBoseOperators()`

discussed in Sec. Defining local operators both create the Fermi phase operator appropriate for this use automatically.

`product`

rule (`mpoterm.ProductTerm`

). Example:

H.AddMPOTerm('product', 'sz', hparam='h_p', weight=-1.0)generates . The operator used in the product must be Hermitian.

`FiniteFunction`

rule (`mpoterm.FiniteFunc`

). Example:

f = [] for ii in range(6): f.append(1.0 / (ii + 1.0)**3) H.AddMPOTerm('FiniteFunction', ['sz', 'sz'], f=f, hparam='J_z', weight=-1.0)generates , where the summation is over all and pairs separated by at least and at most lattice spacings. The range of the potential is given by the number of elements in the list

`f`

. This rule enforces Hermiticity and cares for Fermi phases as in the bond rule case.

`exponential`

rule (`mpoterm.ExpTerm`

). Example:

H.AddMPOTerm('exponential', ['sz', 'sz'], hparam='J_z', decayparam=a, weight=1.0)generates . This rule also enforces Hermiticity and cares for Fermi phases as in the bond rule case.

`InfiniteFunction`

rule (`mpoterm.InfiniteFunc`

). We point out that this term is always fitted with a set of`mpoterm.ExpTerm`

. Examples:

invrcube = lambda x: 1.0/(x*x*x) H.AddMPOTerm('InfiniteFunction', ['sz', 'sz'], hparam='J_z', weight=1.0, func=invrcube, L=1000, tol=1e-9)generates , where the functional form is valid to a distance

`L`

with a residual of at most`tol`

. Similarly,H.AddMPOTerm('InfiniteFunction', ['sz', 'sz'], hparam='J_z', weight=1.0, x=x, y=f, tol=1e-9)generates , where the functional form is determined by the array of values

`f`

and the array of evaluation points`x`

. The distance of validity`L`

is determined from the array`x`

. This rule applies only to monotonically decreasing functions`func`

or`f`

, and so outside the range of validity`L`

the corrections are also monotonically decreasing. This rule enforces Hermiticity and cares for Fermi phases as in the bond rule case.

`MBString`

rule (`mpoterm.MBStringTerm`

). See functions description for usage.

In all of the above expressions, `H`

is an object of the `MPO`

class.
The strings passed to the argument `hparams`

represent tags for variable
Hamiltonian parameters, which allow for easy batching and parallelization
(Sec. Running a parallel static simulation), and also facilitate dynamical processes
The argument `hparams`

is optional, and defaults to the constant
value 1. The optional argument `weight`

specifies a constant factor that
scales the operator and defaults to the value 1, too.

As an example, a template MPO for the Hamiltonian of the hard-core boson model with interactions,

is provided by

```
import MPSPyLib as mps
# Build operators
Operators = mps.BuildBoseOperators(1)
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('site', 'nbtotal', hparam='mu', weight=-1.0)
H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0)
invrcube = lambda x: 1.0/(x*x*x)
H.AddMPOTerm('InfiniteFunction', ['nbtotal', 'nbtotal'], hparam='U', weight=1.0, func=invrcube, L=1000, tol=1e-9)
```

We stress that no assumptions are made about the `hparam`

parameters `'t'`

,
`'mu'`

, or `'U'`

when the MPO is built. The `hparams`

are simply tags
parameterizing the operators which appear in the Hamiltonian. The actual
numerical values of , , and are specified as
part of the dictionary `parameters`

discussed in Sec. Specifying the parameters of a simulation.
These parameters can also be position-dependent, as demonstrated in the example
files in Examples for the openMPS.

One can print operators using the `MPO.MPO.printMPO()`

method. As an
example, the interacting spinless fermion model

is specified and printed as

```
import MPSPyLib as mps
# Build operators
Operators = mps.BuildFermiOperators()
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['fdagger', 'f'], hparam='t', weight=-1.0, Phase=True)
H.AddMPOTerm('bond', ['nftotal', 'nftotal'], hparam='U', weight=1.0)
H.printMPO()
```

The result of this code is

**-1.0 * t * \sum_i fdaggerP_i f_{i+1} + 1.0 * t * \sum_i fP_i fdaggerP_{i+1} + 1.0 * U * \sum_i nftotal_i nftotal_{i+1}**

Note the use of the *phased* operators `fdaggerP`

and `fP`

according to
`Phase=True`

.

As the `'InfiniteFunction'`

rule constructs MPOs by fitting the decaying
function to a series of weighted exponentials, the result of `printMPO`

will contain the data for the exponentials and not information about the
infinite function proper.

Finally, as complex MPOs with `InfiniteFunction`

interactions can be
expensive to generate, one can save the data in an `MPO`

object instance
using the `MPO.WriteMPO()`

function, and then read it in later using
the `MPO.ReadMPO()`

function. The syntax for these is

```
mps.WriteMPO(H,filename)
H = mps.ReadMPO(filename)
```

where `H`

is an object of the `MPO.MPO`

class and `filename`

is
the name of the file where the MPO is written to/read from.

## Defining and reading observables: `Obs.Observables`

¶

### Defining Observables¶

Observables are specified by instantiating an object of the
`Obs.Observables`

class, for example

```
myObservables = Observables(Operators)
```

One then adds observables using the `Obs.Observables.AddObservable()`

method.
Examples of the use of `Obs.Observables.AddObservable()`

are:

`site`

observables. Example :

myObservables.AddObservable('site', 'nbtotal', 'n')will measure at each site and return the result as a length- vector with the tag

`'n'`

. Recall that is either the system size in the case of a finite MPS or the unit cell length in the case of iMPS as specified in the table in the description of`WriteFiles()`

`corr`

observables. Example:

myObservables.AddObservable('corr', ['nbtotal', 'nbtotal'], 'nn')will measure at each independent pair of sites and return the result with the tag

`'nn'`

. For finite-size MPS simulations the result is returned as an matrix, and for infinite-size MPS simulations the result is returned as a length- vector, where is specified by the`SpecifyCorrelationRange`

method, see Sec. Infinite-size ground state search: iMPS. For observables involving odd numbers of fermionic operators on a given site, there is an optional`Phase`

argument such as existed for the`MPO.MPO.AddMPOTerm()`

method. An example ismyObservables.AddObservable('corr', ['fdagger', 'f'], 'spdm', Phase=True)which measures the fermionic correlation function .

`string`

observables. Example (from`Infiniteheisenbergspinone.py`

):

myObservables.AddObservable('string', ['sz', 'sz', 'phaseString'], 'StringOP')will measure at each independent pair of sites and return the result with the tag

`'StringOP'`

. Note that`phaseString`

has been defined to be only in this particular case.

`MPO`

observables. Example:

DefectDensity = mps.MPO(Operators) DefectDensity.AddMPOTerm('bond', ['sz', 'sz']) myStaticObservables.AddObservable('MPO', DefectDensity, 'Defect Density')will measure and return the result with the tag

`'Defect Density'`

. In this way, essentially arbitrary many-body operators may be measured.

- Reduced density matrix, limited to one or two sites at present, can be measured as well. They can be added once as single item or list. An empty list indicates to report the singular values at all bipartitions. If passed as list, entries have to sorted ascending.

myObservables.AddObservables('DensityMatrix_i', []) myObservables.AddObservables('DensityMatrix_ij', [])

It is possible to measure the mutual information matrix without writing all the information about the reduced density matrices to the result files.

myObservables.AddObservables('MI', None)

The second argument,`None`

, is not referenced, but required by the function :py:func:`Obs.Observables.AddObservables.

- The Schmidt values at the bipartitions are another measurement, which gives a more detailed picture than the bond entropy. They can be added once as single item or list. An empty list indicates to report the singular values at all bipartitions.

myObservables.AddObservables('Lambda', [])

Observables for the different algorithms in OSMPS are specified as part of the
`parameters`

dictionary, see Table in `tools.WriteFiles()`

.

### Reading Observables¶

The observables which correspond to static algorithms such as MPS and eMPS,
see Sec. Algorithms of the OpenMPS, are read using the
`Obs.ReadStaticObservables()`

function taking the `parameters`

dictionary as argument. Here, the observables for MPS are specified by the key
`'MPSObservables'`

in `'parameters'`

and, similarly, the key
`'eMPSObservables'`

in `'parameters'`

specifies the observables for eMPS,
see Table in `tools.WriteFiles()`

. The output from
`Obs.ReadStaticObservables()`

is a list of dictionaries. The number of
elements in the list is equal to the number of individual simulations. An
element of this list contains all of the relevant simulation metadata from
`parameters`

as well as the observables specified in `myObservables`

with
the tags given by the `Obs.Observables.AddObservable()`

method. That is,
`Output['n']`

outputs the observable `'n'`

as tagged by
`Obs.Observables.AddObservable()`

for the dictionary `Output`

which is an
element of `Outputs`

. The tags of the various default information returned
from `Obs.ReadStaticObservables()`

can be found in the function
documentation.

In addition to the observables which are specified using the
`Obs.Observables.AddObservable()`

method, the von Neumann entropy of
entanglement of a given site with the rest of the system,

is measured whenever at least one `'site'`

observable is measured. The result
is stored as an -length vector with the tag `'SiteEntropy'`

.
Additionally, the von Neumann entropy of entanglement of bipartition

is measured whenever at least one `'site'`

observable is measured. For
finite-size MPS simulations, the result is stored as an
-length vector with the tag `'BondEntropy'`

. For
infinite-size MPS simulations, `'BondEntropy'`

is always measured, and is a
scalar indicating the entropy of entanglement of the bipartition between any
two adjacent unit cells.

For finite-size systems, the total energy `'energy'`

and the variance
`'variance'`

are always measured.
The quantity `'converged'`

returns a value of `'T'`

or `'F'`

depending
on whether the variance met the desired user-defined tolerance or not,
respectively. For infinite-sized systems, the energy is replaced by
`'energy_density'`

, and the `'convergence'`

is determined from the
orthogonality fidelity as discussed in Sec. Infinite-size ground state search: iMPS. Finally, in all
cases, the maximum bond dimension of the MPS is returned as
`'bond_dimension'`

.

To extract observables from tMPS simulations, the function
`Obs.ReadDynamicObservables()`

taking the `paramaeters`

dictionary
as argument is used. The observables for tMPS are specified by the key
`'DynamicsObservables'`

in `'parameters'`

. The default observables returned
from `Obs.ReadDynamicObservables()`

are given in its function description.
The number of elements of `Outputs`

corresponds to the number of elements of
`parameters`

, but each element of `Outputs`

is itself a list whose length
is the number of time steps. These elements are distinguished by the key
`'time'`

. Time is measured in units of , where
is the unit of energy set by the magnitudes of the
Hamiltonian parameters.

Outputs corresponding to a fixed set of parameters can be extracted from the
output from `Obs.ReadStaticObservables()`

or
`Obs.ReadDynamicObservables()`

using the function
`Obs.GetObservables()`

taking as arguments the `Outputs`

,
`parameter_names`

and `parameter_values`

. We refer to the returned values
as `SmallOutputs`

. Here, `parameter_names`

is a list of parameter names
which are values specified as keys in the `parameters`

passed to
`Obs.ReadStaticObservables()`

and `parameter_values`

is a list of the
values desired for `parameter_names`

. For example, to find only outputs
corresponding to the ground state, one may call

```
GSOutputs = mps.GetObservables(Outputs, 'state', 0)
```

To find outputs corresponding only to the ground state and those simulations
where the Hamiltonian parameter `'t'`

takes the value `0.9`

, one may call

```
GStOutputs = mps.GetObservables(Outputs, ['state', 't'], [0, 0.9])
```

If no matches are found, then an empty list is returned.

## Specification of dynamics¶

A dynamical process simulated with tMPS, details on the algorithms can be
found starting at Sec. Krylov-based time evolution : tMPS, is
specified as an object of the `Dynamics.QuenchList`

class. An object
of this class is obtained as

```
Quenches = mps.QuenchList(H)
```

where `H`

is an `MPO.MPO`

object specifying the MPO of the
Hamiltonian used to time-evolve the state. Dynamical processes, which we will
call quenches, are added using the `Dynamics.QuenchList.AddQuench()`

method with the syntax

```
Quenches.AddQuench(Hparams, time, deltat, funcs, ConvergenceParameters=None, stepsforoutput=1)
```

Here `Hparams`

is a list of the tags of the Hamiltonian
parameters in `H`

which are to be modified during the quench, `time`

is the
total time of this quench, `deltat`

is the desired numerical time step,
`funcs`

is a list of functions specifying how the Hamiltonian parameters given
in `Hparams`

change in time. `ConvergenceParameters`

is used to specify the
time evolution algorithm used and its convergence parameters for this quench.
The available options are the
`convergence.KrylovConvParam`

, `convergence.TEBDConvParam`

,
`convergence.TDVPConvParam`

. and `convergence.LRKConvParam`

.
Finally, `stepsforoutput`

denotes the number time steps between
measurements of the state. The functions in `funcs`

should take a single
argument which is the time. At present, dynamics always begins from the ground
state, and the Hamiltonian used for time evolution must have the same form as
the Hamiltonian used for dynamics. That is, the two Hamiltonians must have the
same Hamiltonian parameters and these parameters must be of the same numerical
type, e.g. if a parameter is site-dependent for statics it must also be
site-dependent for dynamics and vice versa.

It is useful here to consider a specific example. We consider quenching the parameter in the Bose-Hubbard model

linearly from 10 to 1, and then linearly back to 10 over a time . Using the Hamiltonian specification

```
Operators = mps.BuildBoseOperators(5)
Operators['interaction'] = 0.5 * (np.dot(Operators['nbtotal'], Operators['nbtotal']) - Operators['nbtotal'])
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0)
H.AddMPOTerm('site', 'interaction', hparam='U', weight=1.0)
```

this quench process is generated as

```
# Quench function ramping down
def Ufuncdown(t, tau=tau):
return 10.0 + 2.0 * (1.0 - 10.0) * t / tau
# Quench function ramping back up
def Ufuncup(t, tau=tau):
return 1.0 + 2.0 * (10.0 - 1.0) * (t - 0.5 * tau) / tau
Quenches = mps.QuenchList(H)
Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1), [Ufuncdown], ConvergenceParameters=myKrylovConv)
Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1), [Ufuncup], ConvergenceParameters=myKrylovConv)
```

for a given `tau`

and convergence parameters `myKrylovConv`

. We note two
important items regarding the functions `Ufuncdown`

and `Ufuncup`

which
specify how changes in time. The first is that time is measured from
the end of the last quench and not from the beginning of the current quench.
This accounts for the `(t-0.5*tau)`

in the latter quench function. The second
is that `tau`

is included as an optional argument in the functions, and
assigned the desired value. If we had instead used the functions

```
# Quench function ramping down
def Ufuncdown(t):
return 10.0 + 2.0 * (1.0 - 10.0) * t / tau
# Quench function ramping back up
def Ufuncup(t):
return 1.0 + 2.0 * (10.0 - 1.0) * (t - 0.5 * tau) / tau
```

this can lead to incorrect results when looping over tau due to a technical point of python. Generally, any values which are to be looped over in the function definitions should be specified as optional arguments in this fashion.

## Abelian symmetries¶

A particularly important numerical optimization for MPS algorithms is to
respect global symmetries. Respecting a symmetry corresponds to fixing
the quantum number associated with an operator
which commutes with the Hamiltonian . We will call these
quantum numbers *charges*. OSMPS has the ability to simultaneously
conserve an arbitrary number of global U(1) symmetries, which is the most
common type of symmetry found in lattice models. For the Abelian group U(1),
the total charge of a many-body state corresponds to the sum of the charges
of its constituents. That is to say, the total charge operator
may be written

(1)

where is the charge operator for site , which we
will call the *generator* of the symmetry. To conserve charges in OSMPS, one
specifies the generators and their
associated charges. The generators are specified as a list of the keys of the
generators in the dictionary of Operators `Operators`

, see
Sec. Defining local operators. This list of keys is then the value
corresponding to the key `'Abelian_generators'`

in the dictionary
`parameters`

of simulation data, see Sec. Specifying the parameters of a simulation.
The charges are specified as a list of numerical data with the same length as
the list of generators, and is the value corresponding to the key
`'Abelian_quantum_numbers'`

in `parameters`

. We note that the generators
must be diagonal matrices, and an exception will be thrown if this is not the
case. For non-standard symmetries, the operators may be made diagonal by
performing matrix operations on the values of the dictionary `Operators`

as
discussed in Sec. Defining local operators.

In addition to the generators of the symmetry being diagonal, all of the
operators used must be covariant under the group operation. An exception is
thrown if this is not the case. Simply stated, a covariant operator changes
the quantum numbers of a state in a well-defined way. A more complete
definition of covariant operators is given in
Sec. Symmetry-adapted operators of the Developer's manual. For
example, let us consider a system of bosons with the total particle number
conserved. The destruction operator is covariant, as it changes
the particle number by , but the operator
is not covariant, as it does not transform
particle number eigenstates into particle number eigenstates. All operators
which are built by the convenience functions discussed in Sec.
Defining local operators are covariant under the typical group operations, e.g.
all particle operators are covariant under particle number conservation and
all spin operators are covariant under conservation of magnetization.
Operators which are not covariant and are not needed can be removed with the
`RemoveOperator`

function to avoid the aforementioned exception.

In addition to the U(1) symmetry, we support as well symmetries. One example for such a symmetry is the quantum Ising model with the odd and even sector.

At present, conservation of symmetries is only supported for MPS simulations on finite lattices.

## Convergence criteria¶

Convergence criteria are specified as objects inheriting from the class
`convergence.ConvergenceParameters`

. Each algorithm has its own
class with the parameters relevant for the algorithm being used. These are

`convergence.MPSConvParam`

for finite size MPS statics including ground states and excited states.`convergence.iMPSConvParam`

for infinite size MPS statics.- The four time evolution methods for MPS can be adapted with
`convergence.KrylovConvParam`

,`convergence.TEBDConvParam`

,`convergence.TDVPConvParam`

, and`convergence.LRKConvParam`

.

All classes come with default convergence values and those can be used by simply creating an object of the relevant class, e.g.

```
myConv = mps.MPSConvParam()
```

Custom values for the convergence parameters are specified as optional arguments to the object initialization. For example, to set the value of the bond dimension to be 15, one calls

```
myConv = mps.MPSConvParam(max_bond_dimension=15)
```

The convergence parameters contained in `MPSConvParam`

and
`KrylovConvParam`

are discussed in Sec. Variational ground state search and
ond example for the dynamics in Sec. Krylov-based time evolution : tMPS, respectively.

For MPS and eMPS simulations which use the `convergence.MPSConvParam`

class, one can iteratively refine results from previous iterations by
specifying multiple convergence parameters. There are two class methods to do
so. The first is `convergence.MPSConvParam.AddConvergenceParameters()`

,
which has the same arguments as the class instantiation. This includes another
simulation with the new convergence parameters which uses the state obtained
with the previous set of convergence parameters as a variational ansatz. The
second method is
`convergence.ConvParam.AddModifiedConvergenceParameters()`

with the
argumentes refered to as
`copywhich`

, `modifywhich`

, and `whichparameters`

. Here, `copywhich`

is an integer denoting which of the previous convergence parameters are to be
copied, starting from zero. The arguments `modifywhich`

and
`whichparameters`

are lists specifying which parameters are to be modified
and to what values, respectively. For example, to specify that the state is to
be computed with maximum bond dimensions 15, 20, and 25 consecutively, one
would construct the convergence parameters as

```
myConv = mps.MPSConvParam(max_bond_dimension=15)
myConv.AddModifiedConvergenceParameters(0, ['max_bond_dimension'], [20])
myConv.AddModifiedConvergenceParameters(0, ['max_bond_dimension'], [25])
```

For each set of convergence parameters the state is measured, and these outputs
are distinguished by the key `'convergence_parameter'`

, see the key specified
in Table of `Obs.ReadStaticObservables()`

.
The key `'convergence_parameter'`

is an integer beginning
from 0, with increasing integers denoting convergence parameters added later.
For further examples of this functionality, see the example files in
in Sec. Examples for the openMPS.

## Running simulations¶

### Running a serial static simulation¶

As discussed in Sec. Specifying the parameters of a simulation, a list of dictionaries
`parameters`

specifying MPS simulations is converted to Fortran-readable
data using the `tools.WriteFiles()`

method as

```
MainFiles = mps.WriteFiles(parameters, Operators, H)
```

The output `MainFiles`

is a list of the files to be passed to the Fortran
backend. The simulation is run via

```
mps.runMPS(MainFiles)
```

which runs the simulations serially in the order they are specified in the
list `parameters`

. See Sec. Building and installation for information on the
optional arguments to `tools.runMPS()`

appropriate for the type of
build used.

### Running a parallel static simulation¶

OSMPS supports automatic data-parallelism of simulations using MPI (Message Passing Interface). Parallel simulations can be specified in two ways. A list of simulations to be executed as in the serial case or a template where one or multiple parameters can be plugged in from a list. In the following, we explain both variants:

#### Parallel simulations via a list¶

If you already have a list of simulation, you may just switch from the
function `tools.WriteFiles()`

to `tools.WriteMPSParallelFiles()`

.
You need one additional argument, `comp_info`

, which is a dictionary
containing the setting for a computer cluster with the only required argument
the number of nodes. For a better tuning more parameters are available. The
example is then

```
comp_info = {'nodes' : '4'}
MainFiles = mps.WriteFiles(parameters, Operators, H, comp_info)
```

This will return the name of script to be submitted to the computer cluster.
This is step replaces the previous call to `tools.runMPS()`

. In case
you want to run the job with MPI on a local machine, you will inside the
script a line containing **Execute_MPSParallelMain**. This line can
be copied to you command line and starts the MPI simulations.

#### Parallel simulations via a template¶

We assume that there exist some set of parameters which one wishes to parallelize over. For example, one may wish to parallelize the XXZ model in a magnetic field

over the anisotropy and the magnetic field . One begins
by building a template dictionary `parameters_template`

for the simulations
which includes all information except for the specification for
and as discussed in Sec. Specifying the parameters of a simulation. The specification of
the parameters and comes through two lists. The first,
called `iterparams`

, is a list of the tags of the Hamiltonian parameters to
be parallelized over. The second, called `iters`

, is a list whose elements
are a list of the numerical values for the corresponding parameters in
`iterparams`

. For our example of the XXZ model, these two vectors might be
specified as

```
Deltamin = 0.1
Deltamax = 2.0
Deltavals = 50
Deltaiter = np.linspace(Deltamin, Deltamax, Deltavals)
hmin = 0.0
hmax = 1.0
hvals = 20
hiter = np.linspace(hmin, hmax, hvals)
iterparams = ['Delta', 'h']
iters = [Deltaiter, hiter]
```

This would parallelize over the 50 values of linearly spaced
from 0.1 to 2 and the 20 values of linearly spaced from 0 to 1 for
a total of 1000 points. Note that the tags we use here in `iterparams`

should
match those we use in building the Hamiltonian MPO, see Sec. Defining the Hamiltonian and other many-body operators.
These lists, along with the template dictionary of other metadata, are passed
to the function `Paralleltools.WriteMPSParallelTemplate()`

, with syntax

```
parameters = mps.WriteMPSParallelTemplate(parameters_template, Operators, H, comp_info, iterparams, iters)
```

Here, `comp_info`

is a dictionary specifying metadata related to the parallel
environment. The function `Paralleltools.WriteMPSParallelTemplate()`

writes the files required by the Fortran backend for each individual simulation
as well as two files `*jobpool.dat`

and `*progress`

, where `*`

is the
`'job\_ID'`

specified in the `parameters_template`

dictionary. The jobpool
file contains the main files which specify the individual simulations to be run
together with unique integer tags, and the progress file indicates which of
these tasks has already been completed. Finally,
`Paralleltools.WriteMPSParallelTemplate()`

writes a pbs or slurm script
which will submit a parallel job to your computer's queueing system. If a job
does not finish all of the tasks in the jobpool, calling
`Paralleltools.WriteMPSParallelTemplate()`

will create a new job pool
containing only simulations which have not been run previously.

An alternate functionality of the parallel code is to pass in a list of
parameters to `parameters_template`

.

`Paralleltools.WriteMPSParallelTemplate()`

will automatically parallelize
over all of the simulations in ```
parameters_template`. A list of Hamiltonians
may also be passed to ``H
```

, but the length of `H`

and
`parameters_template`

must be the same. Finally, the two functionalities
may be combined by specifying both a list of `parameters_template`

as well
as parameters to iterate over via `iterparams`

and `iters`

.

## Putting it all together¶

In this section we will analyze the simplest example file given in the
`Examples`

subdirectory of OSMPS, `SpinlessFermions.py`

. The file begins
by importing *MPSPyLib*, *NumPy*, and the module *sys*

1 2 3 | ```
import MPSPyLib as mps
import numpy as np
import sys
``` |

Inside the function `main`

we start building the spinless operators

1 2 | ```
# Build operators
Operators = mps.BuildFermiOperators()
``` |

Now, `Operators`

contains the fermion operators given in the table of
`Ops.BuildFermiOperators()`

. We now build the Hamiltonian which consists
only of tunneling,

(2)

This is done as

1 2 3 | ```
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['fdagger', 'f'], hparam='t', weight=-1.0, Phase=True)
``` |

The specification `'bond'`

denotes the sum over nearest-neighbors in
Eq. (2) and `Phase=True`

denotes that a Fermi phase should be
included, as this term involves an odd number of fermion operators at each
site, see Sec. Defining the Hamiltonian and other many-body operators. We next define the observables we wish to measure

1 2 3 4 5 | ```
# Observables
myObservables = mps.Observables(Operators)
# Correlation functions
myObservables.AddObservable('corr', ['fdagger','f'], 'bspdm')
myObservables.AddObservable('corr', ['fdagger','f'], 'spdm', Phase=True)
``` |

Here, the observable we have called `'spdm'`

is the fermionic
single-particle density matrix, , and the observable we have called `'bspdm'`

is a single-particle
density matrix in which the operators anticommute on a
given site but commute when not on the same site. The difference is in
including `Phase=True`

in the observable specification. We will see the
differences in the two correlation functions when we run the simulation. Next,
we specify the convergence criteria for our simulation as

1 2 | ```
# Convergence data
myConv = mps.MPSConvParam(max_bond_dimension=30, max_num_sweeps=2)
``` |

Now, we specify the simulation parameters as

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | ```
t = 1.0
L = 10
N = 5
# Define statics
parameters = [{
# Directories
'simtype' : 'Finite',
'job_ID' : 'SpinlessFermions_',
'unique_ID' : 'L_'+str(L)+'N'+str(N),
'Write_Directory' : 'TMP/',
'Output_Directory' : 'OUTPUTS/',
# System size and Hamiltonian parameters
'L' : L,
't' : t,
'verbose' : 2,
'logfile' : True,
# Specification of symmetries and good quantum numbers
'Abelian_generators' : ['nftotal'],
'Abelian_quantum_numbers' : [N],
'MPSObservables' : myObservables,
``` |

`L`

is the number of lattice sites and `t`

is the tunneling energy. Note
that it is only in the dictionary `parameters`

that we specify a value for
`t`

; we do not need to give an actual value for `t`

when we are building
the Hamiltonian. The lines

1 2 | ```
'logfile' : True,
# Specification of symmetries and good quantum numbers
``` |

specify that the total number of fermions is conserved, and is set to be `N`

.
We write Fortran-readable files and run them with the lines if we are not in
the post processing mode:

1 2 3 4 5 6 7 | ```
}]
# Write Fortran-readable main files
MainFiles = mps.WriteFiles(parameters, Operators, H,
PostProcess=PostProcess)
if(not PostProcess):
``` |

The observables are read in using

1 2 3 4 | ```
mps.runMPS(MainFiles, RunDir='./')
# Postprocessing
# --------------
``` |

Because we have only a single simulation in this case, corresponding to the
ground state with a single set of convergence parameters, the relevant outputs
are contained in `Outputs[0]`

, the first element of the list `Outputs`

, see
Sec. Defining and reading observables: Obs.Observables. The matrix representing the
single-particle density matrix is `Outputs[0]['spdm']`

, as we specified the
tag `'spdm'`

for this observable in line 14. We now diagonalize the
single-particle density matrix and its *phaseless* counterpart `'bspdm'`

to
find their eigenvalues using the numpy function `eigh`

and print the
eigenvalues to the screen with

1 2 3 4 5 6 7 8 9 | ```
Outputs = mps.ReadStaticObservables(parameters)
spdm = Outputs[0]['spdm']
spdmeigs, U = np.linalg.eigh(spdm)
bspdm = Outputs[0]['bspdm']
bspdmeigs, U = np.linalg.eigh(bspdm)
print('Eigenvalues of <f^{\dagger}_i f_j> with Fermi phases', spdmeigs)
print('Eigenvalues of <f^{\dagger}_i f_j> without Fermi phases', bspdmeigs)
``` |

The return state indicates the end of the main function. The last lines call
our main function in the case the file is call from command line and handle
the post processing argument from the command line. If it was included as
module the `main`

function would not be executed.

1 2 3 4 5 6 7 8 9 | ```
if(__name__ == '__main__'):
# Check for command line arguments
Post = False
for arg in sys.argv[1:]:
key, val = arg.split('=')
if(key == '--PostProcess'): Post = (val == 'T') or (val == 'True')
``` |

Calling

```
osmps@manual:~$ python SpinlessFermions.py
```

in the terminal gives the output

```
./Execute_MPSMain TMP/SpinlessFermions_L_10N5Main.nml
Eigenvalues of <f^{\dagger}_i f_j> with Fermi phases [ -4.56780243e-16 -4.34812894e-17 2.50202620e-16
5.99741719e-13 6.00558567e-13 1.00000000e+00 1.00000000e+00 1.00000000e+00
1.00000000e+00 1.00000000e+00]
Eigenvalues of <f^{\dagger}_i f_j> without Fermi phases [ 0.05067985 0.06159355 0.09456777 0.12959356
0.17094556 0.35616646 0.45241313 0.60737271 0.91929511 2.1573723 ]
```

We see that the properly phased single-particle density matrix produces a sharp Fermi distribution, as expected. Contrariwise, the single-particle density matrix without the Fermi phases displays a macroscopically large eigenvalue.

Footnotes