.. _sec_simulationdetails:

====================
 Simulation Details
====================

.. \label{sec:Include}
.. _sec_including_libraries:

Including the required Python libraries
=======================================

Python libraries are included in a file by placing

.. code-block:: python

   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 <http://docs.python.org/tutorial/classes.html#python-scopes-and-namespaces>`_
where the contents of this module are placed. As an example, if we wanted to
include MPSPyLib and call the function :py:func:`Ops.BuildSpinOperators` from
this library, we would use the code

.. code-block:: python

   import MPSPyLib as mps
   Operators = mps.BuildSpinOperators()

Note the use of the namespace identifier ``mps`` to specify that
:py:func:`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.


.. _sec_parameters:

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 <http://www.tutorialspoint.com/python/python_dictionary.htm>`_, 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

.. code-block:: python

   mydict = {'name' : 'fred'}

A dictionary is specified by curly braces, :math:`\{ \ldots \}`, 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, :math:`\left[ \ldots \right]`, 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 :py:func:`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. :ref:`sec_mpo`. The simulation specification ``parameters``
is converted into simulation data which is readable by the Fortran backend using
the function

.. code-block:: python

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

where ``Operators`` specifies the operators defining the on-site Hilbert space,
see Sec. :ref:`sec_deflocalops`, and ``H`` specifies the Hamiltonian, see
Sec. :ref:`sec_mpo`. 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.


.. \label{sec:Opsdpy}
.. _sec_deflocalops:

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.
:ref:`sec_abelian`). This means that the set of operators is stored as a
dictionary, call it ``Operators``, and

.. code-block:: python

   Operators['tag']

returns the matrix representation of the operator ``'tag'`` as a 
`NumPy array <http://docs.scipy.org/doc/numpy/reference/generated/numpy.array.html>`_.
The default operators built by the convenience functions are listed in a
table inside of the function description. The options are spin operators
from :py:func:`Ops.BuildSpinOperators`, bosonic operators defined in
:py:func:`Ops.BuildBoseOperators`, fermionic operators created by the
function :py:func:`Ops.BuildFermiOperators`, or operators for a Bose-Fermi
system :py:func:`Ops.BuildFermiBoseOperators` building all operators present
in the bosonic and fermionic set. The total number operator :math:`\hat{n}`
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.

.. code-block:: python

   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

.. math::

   \hat{U} &= \frac{1}{2} \left( \hat{n}^2 - \hat{n} \right) \, ,

and call it ``'Interaction'``. This may be achieved as

.. code-block:: python

   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
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html>`_.

.. Operators can be removed using the  ``RemoveOperator`` convenience
   function, though this should rarely be necessary. (It is implemented
   but still might cause trouble ...)


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 :py:func:`Ops.BuildBoseOperators` or
:py:func:`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 :math:`2s`, with :math:`s` 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 :py:func:`Ops.BuildBoseOperators` or
:py:func:`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 :math:`2s`, with :math:`s` the spin. In
:py:func:`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 :py:func:`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

.. code-block:: python

   import MPSPyLib as mps
   import numpy as np

   Operators = mps.BuildBoseOperators(Nmax=2)
   print(Operators['b'])

Likewise, we can print all the operators as

.. code-block:: python

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


.. _sec_mpo:

Defining the Hamiltonian and other many-body operators
======================================================

.. automodule:: mpo


.. _sec_observables:

Defining and reading observables: :py:class:`obsterms.Observables`
==================================================================

.. automodule:: obsterms

Reading Observables
-------------------

The observables which correspond to static algorithms such as MPS and eMPS,
see Sec. :ref:`sec_algorithms`, are read using the
:py:func:`obsterms.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  :py:func:`tools.WriteFiles`.  The output from
:py:func:`obsterms.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 :py:func:`obsterms.Observables.AddObservable` method. That is,
``Output['n']`` outputs the observable ``'n'`` as tagged by
:py:func:`obsterms.Observables.AddObservable` for the dictionary ``Output`` which is an
element of ``Outputs``. The tags of the various default information returned
from :py:func:`obsterms.ReadStaticObservables` can be found in the function
documentation.
    
In addition to the observables which are specified using the
:py:func:`obsterms.Observables.AddObservable` method, the von Neumann entropy of
entanglement of a given site :math:`i` with the rest of the system,

.. math::

   S_{\mathrm{site},i}&\equiv -\mathrm{Tr}\left[\hat{\rho}_i\log\hat{\rho}_i\right]\, ,\;\;\; \hat{\rho}_i\equiv \mathrm{Tr}_{j\ne i}|\psi\rangle\langle \psi|

is measured whenever not deactivated. The result
is stored as an :math:`L`-length vector with the tag ``'SiteEntropy'``.
Additionally, the von Neumann entropy of entanglement of bipartition :math:`i`

.. math::

   S_{\mathrm{bond},i}&\equiv -\mathrm{Tr}\left[\hat{\rho}_i\log\hat{\rho}_i\right]\, ,\;\;\; \hat{\rho}_i\equiv \mathrm{Tr}_{j\ge i}|\psi\rangle\langle \psi|

is measured whenever not deactivated. For
finite-size MPS simulations, the result is stored as an
:math:`\left(L+1\right)`-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'`` :math:`\equiv\langle \hat{H}^2-E^2\rangle` 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. :ref:`sec_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
:py:func:`obsterms.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 :py:func:`obsterms.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 :math:`\hbar/E_{\mathrm{unit}}`, where
:math:`E_{\mathrm{unit}}` is the unit of energy set by the magnitudes of the
Hamiltonian parameters.


.. _sec_dynamics:

Specification of dynamics
=========================

A dynamical process simulated with tMPS, details on the algorithms can be
found starting at Sec. :ref:`sec_krylov`, is
specified as an object of the :py:class:`Dynamics.QuenchList` class. An object
of this class is obtained as

.. code-block:: python

   Quenches = mps.QuenchList(H)

where ``H`` is an :py:class:`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 :py:func:`Dynamics.QuenchList.AddQuench`
method with the syntax

.. code-block:: python
		    
   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
:py:class:`convergence.KrylovConvParam`, :py:class:`convergence.TEBDConvParam`,
:py:class:`convergence.TDVPConvParam`. and :py:class:`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 :math:`U` in the Bose-Hubbard model

.. math::

   \hat{H}&=-t\sum_{\langle i,j\rangle}\left[\hat{b}_i^{\dagger}\hat{b}_j+\mathrm{h.c.}\right]+U\sum_{i}\frac{\hat{n}_i\left(\hat{n}_i-1\right)}{2}

linearly from 10 to 1, and then linearly back to 10 over a time :math:`\tau`.
Using the Hamiltonian specification

.. code-block:: python

   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

.. code-block:: python

   # 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 :math:`U` 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

.. code-block:: python

   # 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 <http://stackoverflow.com/questions/233673/lexical-closures-in-python>`_.
Generally, any values which are to be looped over in the function definitions
should be specified as optional arguments in this fashion.


.. _sec_abelian:

Abelian symmetries
==================

A particularly important numerical optimization for MPS algorithms is to
respect global symmetries.  Respecting a symmetry corresponds to fixing
the quantum number :math:`q` associated with an operator :math:`\hat{Q}`
which commutes with the Hamiltonian :math:`\hat{H}`. 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 :math:`\hat{Q}`
may be written


.. math::
   :label: eqgens

   \hat{Q} &= \sum_{i} \hat{Q}_i \, ,


where :math:`\hat{Q}_i` is the charge operator for site :math:`i`, which we
will call the *generator* of the symmetry. To conserve charges in OSMPS, one
specifies the generators :math:`\left\{\hat{Q}_i\right\}` 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. :ref:`sec_deflocalops`. This list of keys is then the value
corresponding to the key ``'Abelian_generators'`` in the dictionary
``parameters`` of simulation data, see Sec. :ref:`sec_parameters`.
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. :ref:`sec_deflocalops`.

.. uncertain reference eq:SAOp for covariant operators, dlarue 20140228
   I've changed the \ref{eq:SAOp} reference to the sec:SymmetryAdaptedOperators,
   which I believe is the intended reference. The pdf should reference
   Section 12.2., aglick 20140301

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. :ref:`sec_symmetryadpatedoperators` of the Developer's manual. For
example, let us consider a system of bosons with the total particle number
conserved. The destruction operator :math:`\hat{b}` is covariant, as it changes
the particle number by :math:`-1`, but the operator
:math:`\hat{b}+\hat{b}^{\dagger}` 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.
:ref:`sec_deflocalops` 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 :math:`\mathbb{Z}_{2}`
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.


.. _sec_conv:

Convergence criteria
====================

Convergence criteria are specified as objects inheriting from the class
:py:class:`convergence.ConvergenceParameters`. Each algorithm has its own
class with the parameters relevant for the algorithm being used. These are
    
* :py:class:`convergence.MPSConvParam` for finite size MPS statics including
  ground states and excited states.
* :py:class:`convergence.iMPSConvParam` for infinite size MPS statics.
* The four time evolution methods for MPS can be adapted with
  :py:class:`convergence.KrylovConvParam`,
  :py:class:`convergence.TEBDConvParam`,
  :py:class:`convergence.TDVPConvParam`,
  and :py:class:`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.  

.. code-block:: python

   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

.. code-block:: python

   myConv = mps.MPSConvParam(max_bond_dimension=15)

The convergence parameters contained in :py:class:`MPSConvParam` and
:py:class:`KrylovConvParam` are discussed in Sec. :ref:`sec_vmps` and
ond example for the dynamics in Sec. :ref:`sec_krylov`, respectively.

For MPS and eMPS simulations which use the :py:class:`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 :py:func:`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
:py:func:`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

.. code-block:: python

   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 :py:func:`obsterms.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. :ref:`sec_examples`.


.. _sec_running:

Running simulations
===================

Running a serial static simulation
----------------------------------

As discussed in Sec. :ref:`sec_parameters`, a list of dictionaries
``parameters`` specifying MPS simulations is converted to Fortran-readable
data using the :py:func:`tools.WriteFiles` method as

.. code-block:: python

   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

.. code-block:: python

   mps.runMPS(MainFiles)

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


.. _sec_parallel:

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 :py:func:`tools.WriteFiles` to :py:func:`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

.. code-block:: python

   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 :py:func:`tools.runMPS`. In case
you want to run the job with MPI on a local machine, you will inside the
script a line containing :command:`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

.. math::

   \hat{H} = \sum_{\langle i,j \rangle} \left(\hat{S}_i^ + \hat{S}_j^- + \mathrm{h.c.} + \Delta \hat{S}^z_i \hat{S}^z_j \right)- h \sum_i \hat{S}^z_i

over the anisotropy :math:`\Delta` and the magnetic field :math:`h`. One begins
by building a template dictionary ``parameters_template`` for the simulations
which includes all information except for the specification for :math:`\Delta`
and :math:`h` as discussed in Sec. :ref:`sec_parameters`. The specification of
the parameters :math:`\Delta` and :math:`h` 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

.. code-block:: python

   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 :math:`\Delta` linearly spaced
from 0.1 to 2 and the 20 values of :math:`h` 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. :ref:`sec_mpo`.
These lists, along with the template dictionary of other metadata, are passed
to the function :py:func:`Paralleltools.WriteMPSParallelTemplate`, with syntax

.. code-block:: python

   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 :py:func:`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,
:py:func:`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
:py:func:`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``.

:py:func:`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``.


.. _sec_example_spinlessfermions:

Putting it all together
=======================

In this section we will analyze the simplest example file given in the
``Examples`` subdirectory of OSMPS, ``04_SpinlessFermions.py``. The file begins
by importing *MPSPyLib*, *NumPy*, and the module *sys*

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 1-3
   :linenos:

Inside the function ``main`` we start building the spinless operators

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 18-19
   :linenos:

Now, ``Operators`` contains the fermion operators given in the table of
:py:func:`Ops.BuildFermiOperators`. We now build the Hamiltonian which consists
only of tunneling,

.. math::
   :label: eqspf

   \hat{H}&=-t\sum_{\langle i,j\rangle}\left[\hat{f}_i^{\dagger}\hat{f}_j+\mathrm{h.c.}\right]\, .

This is done as

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 21-23
   :linenos:

The specification ``'bond'`` denotes the sum over nearest-neighbors in
Eq. :eq:`eqspf` 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. :ref:`sec_mpo`. We next define the observables we wish to measure

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 25-29
   :linenos:

Here, the observable we have called ``'spdm'`` is the fermionic
single-particle density matrix, :math:`\langle \hat{f}_i^{\dagger}\hat{f}_j
\rangle`, and the observable we have called ``'bspdm'`` is a single-particle
density matrix :math:`\langle \hat{\tilde{f}}_i^{\dagger}\hat{\tilde{f}}_j
\rangle` in which the operators :math:`\hat{\tilde{f}}_i` 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

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 31-32
   :linenos:

Now, we specify the simulation parameters as

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 34-54
   :linenos:

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

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 50-51
   :linenos:

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:

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 56-62
   :linenos:

The observables are read in using

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 64-67
   :linenos:

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. :ref:`sec_observables`. The :math:`L\times L` 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

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 69-77
   :linenos:

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.

.. literalinclude:: ../../../Examples/04_SpinlessFermions.py
   :lines: 80-88
   :linenos:

Calling

.. code-block:: bash

   osmps@manual:~$ python 04_SpinlessFermions.py

in the terminal gives the output

.. code-block:: bash

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

.. rubric:: Footnotes

.. Note that the phased Fermi operators ``fP`` and ``fdaggerP`` are
       only constructed if they are linearly independent of ``f`` and
       ``fdagger``, i.e. not for hard-core fermions.