.. _Examples:

########
Examples
########

We illustrate the usage using various examples that include the different :ref:`Operator types<Operator>` (Lindbladian/Hamiltonian with (**U1**) and without (**NoU1**) number conserving) as well as different problems such as time evaluation and ground state searches.


All examples are attached in the repository and provide a great starting point for **your project**.
You can create a build folder in the respective directories and call :ref:`CMake<Config_and_CMake>` to build the source files:

.. code-block:: cpp

   cd your/example
   mkdir build & cd build
   cmake .. & make
   ./main -some_options ...

Overview
--------

We provide multiple examples handling the different libraries that can be used with our implementation as well as different problems such as ground state (Hamiltonian) and steady states (Lindbladian) searches, and imaginary and real time evolutions.
The following list provides a brief overview about all available examples:

* :ref:`Examples_Basis`: Testing out **DanceQ** without extension

    * A simple script that illustrate the usage of :ref:`State class<State>` and :ref:`BasisU1 class<BasisU1>`

* :ref:`Examples_Operator_class`: Testing out **DanceQ** without extension

    * A simple script that illustrate usage of the :ref:`Operator class<Operator>`

* :ref:`Examples_Eigen`: Full diagonalization with :ref:`Eigen<Eigen_setup>`

    * This example computes **all** eigenvalues of a Lindbladian with number conservation
    
    * It is based on :ref:`Eigen<Eigen_setup>` and can be linked to `MKL <https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html>`_

    * Code is compatible :ref:`openMP<openMP_docs>`

* :ref:`Examples_ShellMatrix`: Ground state energy computation with our own :ref:`Lanczos<Functions_lanczos>` implementation

    * This example uses the :ref:`ShellMatrix<ShellMatrix>` for matrix-free multiplication to compute the ground state of Hamiltonian with number conservation
    
    * We provide :ref:`Lanczos<Functions_lanczos>` implementation that only stores **two states** to reduce the memory consumption

    * Code is compatible :ref:`openMP<openMP_docs>` and :ref:`MPI<MPI_docs>`

* :ref:`Examples_SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution`: **Imaginary** time evolution with :ref:`Petsc and Slepc<Petsc_setup>`

    * This example computes the **imaginary** time evolution :math:`e^{-\beta H}\vert\Psi\rangle` for a Hamiltonian system without number conservation

    * The code is based on :ref:`Petsc and Slepc<Petsc_setup>` that distributes the sparse operator across multiple ranks

    * Code is compatible :ref:`MPI<MPI_docs>`


* :ref:`Examples_SparseMatrix_petsc_Hamiltonian_real_time_evolution`: **Real** time evolution with :ref:`Petsc and Slepc<Petsc_setup>`

    * This example computes the **real** time evolution :math:`e^{-it H}\vert\Psi\rangle` for a Hamiltonian system with number conservation

    * The code is based on :ref:`Petsc and Slepc<Petsc_setup>` that distributes the sparse operator across multiple ranks

    * Code is compatible :ref:`MPI<MPI_docs>`

* :ref:`Examples_SparseMatrix_petsc_Lindbladian`: Open-system dynamics with :ref:`Petsc and Slepc<Petsc_setup>`

    * This example computes the real time evolution :math:`e^{tH}\vert\rho\rangle` for a Lindbladian system with number conservation

    * The code is based on :ref:`Petsc and Slepc<Petsc_setup>` that distributes the sparse operator across multiple ranks

    * Code is compatible :ref:`MPI<MPI_docs>`

* :ref:`Examples_ShellMatrix_petsc_Hamiltonian`: Ground state computation with :ref:`Petsc and Slepc<Petsc_setup>`

    * This example computes the ground state of a Hamiltonian system with number conservation

    * The code is based on :ref:`Petsc and Slepc<Petsc_setup>` that distributes the sparse operator across multiple ranks

    * Code is compatible :ref:`MPI<MPI_docs>`


* :ref:`Examples_ShellMatrix_petsc_Lindbladian`: Steady-state computation with :ref:`Petsc and Slepc<Petsc_setup>`

    * This example computes the steady states and slowest decay rates of a Lindbladian system without number conservation and non-hermitian jump operators

    * The code is based on :ref:`Petsc and Slepc<Petsc_setup>` that distributes the sparse operator across multiple ranks

    * Code is compatible :ref:`MPI<MPI_docs>`


.. _Examples_Basis:

State and Basis class
---------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/Basis.rst


The first example introduces the :ref:`State<State>` and :ref:`Basis class<Basis_General>`.
The core of :ref:`DanceQ<DanceQ>` is the :ref:`State class<State>` that **bitwise** encodes a product state and provides a great functionality.
Most functions act on an interval defined by **start** and **end**.

.. code-block:: cpp

    /* Default constructor, all sites are set to zero */
    State state;

    /* Sets the last sites to 1 */
    state.set_range(1, /*start=*/0, /*end=*/L);

    /* Sets the last state to Q-1 */
    state.set_state(0,L-1);

    /* Output of the state */
    std::cout << "[main] - state: " << state.get_string(/*start=*/0UL, /*end=*/L) << std::endl;

    /* Output of the first site */
    std::cout << "[main] - first site: " << state[0] << std::endl;

    /* Output of the last site */
    std::cout << "[main] - last site: " << state[L-1] << std::endl;


Further, the example illustrates the use of the :ref:`Basis class<Basis_General>` with the :ref:`Basis Iterator<Iterator>`.
The :ref:`Basis class<BasisU1>` is built on top of the :ref:`Container class<Container>` that built on top of the :ref:`State class<State>`.
It is constructed by defining the system size **L**, the particle sector **n**, and the number of subsystems **Npart**.
The :ref:`Iterator class<Iterator>`, **it** below, allows an easy handling of the :ref:`Basis<Basis_General>`:

.. code-block:: cpp

    /* Construct the Basis of size L, with total particle number n, and Npart subsystems */
    Basis basis(L,n,Npart);

    /* Prints Basis information */
    basis.info();

    /* Simple iteration using the Basis iterator */
    for(auto it=basis.begin(); it<basis.end(); ++it){

        /* Retrieving the index from the iterator */
        uint64_t current_index = it.get_index();

        /* Output of the state; it-> points to the State class */
        std::cout << "[main] - |Psi_" << current_index << "> = " << it->get_string(/*start=*/0UL, /*end=*/L) << std::endl;

    }


.. _Examples_Operator_class:

Operator class
--------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/Operator.rst


The :ref:`Operater class<Operator>` can be either built using :ref:`Basis class<BasisU1>` or the :ref:`State class<State>`.
In the latter case, the :ref:`operator<Operator>` acts on the full Hilbert space with dimension :math:`\mathrm{Q}^L` without number conservation:

.. code-block:: cpp

    /* Definition of the State class */
    using State=danceq::internal::State<MaxSites,Q>;

    /* Definition of the Basis class */
    using Basis=danceq::internal::BasisU1<State>;

    /* Definition of the Operator class from the Basis class */
    using MyOperator_U1=danceq::internal::Operator<Basis,double>;

    /* Definition of the Operator class from the State class */
    using MyOperator_NoU1=danceq::internal::Operator<State,double>;


:ref:`DanceQ<DanceQ>` provides a user-friendly interface to define the Hamiltonian.
This example defines an one-dimensional XXX-Heisenberg model with periodic boundary conditions:

.. math::

   H = \sum_{i=0}^{L-1} \frac{1}{2}\left(S^+_iS^-_{i+1} +S^-_iS^+_{i+1}\right) + S^z_iS^z_{i+1}


.. code-block:: cpp

    MyOperator_U1 H(L,n);

    for(uint64_t i = 0; i < L; i++){
        H.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
        H.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
        H.add_operator(1., {i,(i+1)%L}, {"Sz","Sz"});
    }

The Hamiltonian can be retrieved in a dense format using a two-dimensional ``std::vector`` with :ref:`create_DenseMatrix()<Operator_Functions_Frontends_create_DenseMatrix>` or in a sparse format, see :ref:`SparseMatrix<SparseMatrix>`, with :ref:`create_SparseMatrix()<Operator_Functions_Frontends_create_SparseMatrix>`.

.. code-block:: cpp

    auto H_dense = H.create_DenseMatrix();
    auto H_sparse = H.create_SparseMatrix();

Both functions are compatible with :ref:`openMP<openMP_docs>`.
This can be activated with ``cmake -DOMP=ON ..``.

.. _Examples_Eigen:

Eigen full spectrum
-------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/Eigen.rst


Functions in this :ref:`example<Examples_Eigen_Code>` are available if the ``EIGEN_DIR`` is set correctly.
A detailed guide to set up `Eigen <https://eigen.tuxfamily.org/index.php?title=Main_Page>`_ can be found :ref:`here<Eigen_setup>`.
:ref:`Eigen<Eigen_setup>` should be used for dense computations such as determining the full spectrum of an :ref:`operator<Operator>`.
In this example, we look at an one-dimensional XXX-Heisenberg model with periodic boundary conditions and an alternating field in :math:`S^z_i`.

.. math::

   H = \sum_{i=0}^{L-1} \frac{1}{2}\left(S^+_iS^-_{i+1} +S^-_iS^+_{i+1}\right) + S^z_iS^z_{i+1} + (-1)^i S^z_i

The Lindbladian further includes a jump operator acting on the first site:

.. math::

   L[\rho] = -i[H,\rho] + S_0^z\rho {S_0^z}^\dagger - \frac{1}{2}\left\{{S_0^z}^\dagger S_0^z,\rho\right\}


.. code-block:: cpp

	/* Operator class */
	Lindbladian_U1 Lind(L,n);

	uint64_t dim = Lind.get_dim();
	std::cout << "[main] - dim = " << dim << std::endl;

	/* Simply XXX-Heisenberg for the unitary dynamics  */
	for(uint64_t i = 0; i < L; i++){
		Lind.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
		Lind.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
		Lind.add_operator(1., {i,(i+1)%L}, {"Sz","Sz"});
		Lind.add_operator(std::pow(-1.,i), {i}, {"Sz"});
	}

	/* Jump operator acting on the first site */
	Lind.add_jump_operator(1., {0}, {"Sz"});

The :ref:`Eigen matrix<Eigen_setup>` can be retrieved with :ref:`create_EigenDenseMatrix()<Operator_Functions_Frontends_create_EigenDenseMatrix>`.
The functions is compatible with :ref:`openMP<openMP_docs>`.

.. code-block:: cpp

	auto Lind_matrix = Lind.create_EigenDenseMatrix();

The returned matrix can be used in all `Eigen routines <https://eigen.tuxfamily.org/dox-devel/classEigen_1_1SelfAdjointEigenSolver.html>`_.
The following code determines all eigen values:

.. code-block:: cpp

    /* Non-hermitian complex solver */
    Eigen::ComplexEigenSolver<decltype(Lind_matrix)> solver;

    /* Computation without eigenvectors */
    solver.compute(Lind_matrix, false);

    auto eigenvalues = solver.eigenvalues();	

If :ref:`MKL<Eigen_setup>` is enabled with ``cmake -DMKL=ON -DOMP=ON ..``, the code can parallelized (4 cores in this case) with:

.. code-block:: console

    export OMP_NUM_THREADS=4
    export MKL_NUM_THREADS=4

    ./main -L 10 -n 5

.. _Examples_ShellMatrix:

ShellMatrix ground state
-------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/ShellMatrix.rst


We determine the ground state energy using our own **matrix-free** implementation given in the :ref:`ShellMatrix class<ShellMatrix>`.
This is a native **cpp** implementation supporting distributed (:ref:`MPI<MPI_docs>`) and shared (:ref:`openMP<openMP_docs>`) memory.
:ref:`CMake<Config_and_CMake>` allows you to simply switch between both by specifying ``cmake -DMPI=ON -DOMP=OFF ..`` or ``cmake -DMPI=OFF -DOMP=ON ..``.




Given an :ref:`Operator<Operator>`, the :ref:`ShellMatrix<ShellMatrix>` provides a **matrix-vector** multiplication (**without storing the matrix**) that can be used to generate `Krylov subspaces <https://en.wikipedia.org/wiki/Krylov_subspace>`_ to compute ground states and other properties.
States are represents by the :ref:`Vector class<Vector>` that -- as well -- supports :ref:`MPI<MPI_docs>` and :ref:`openMP<openMP_docs>`.


In the example, we focus on an one-dimensional XXZ-Heisenberg model with open boundary conditions and a tilted field:

.. math::

   H = \sum_{i=0}^{L-2} \frac{1}{2}\left(S^+_iS^-_{i+1} +S^-_iS^+_{i+1}\right) + J_z S^z_iS^z_{i+1} + \sum_{i=0}^{L-1} i S^z_i


.. code-block:: cpp

    Hamiltonian_U1 H(L,n);
    
    for(uint64_t i = 0; i < L-1; i++){
        H.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
        H.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
        H.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
        H.add_operator(static_cast<double>(i)., {i}, {"Sz"});
    }

We further set up an observable that measures the magnetization on the first site.
The :ref:`ShellMatrix<ShellMatrix>` is simply generated with :ref:`create_ShellMatrix()<Operator_Functions_Frontends_create_ShellMatrix>`:

.. code-block:: cpp

	Hamiltonian_U1 O(L,n); 

	O.add_operator(1., {0}, {"Sz"});

	auto O_matrix = O.create_ShellMatrix();



We have implemented the `Lanczos algorithm <https://en.wikipedia.org/wiki/Lanczos_algorithm>`_ to compute the ground state and excitations with :ref:`ShellMatrix class<ShellMatrix>`.
It is simply called by passing the Hamiltonian to the function :ref:`lanczos(H,...)<Functions_lanczos>`.

.. code-block:: cpp

    /* Vector to collect the eigen states */
    std::vector<Vector> states;

    /* Returns the eigen values */
    auto data = lanczos(H, number_of_states, &states /* returns eigen states if this is not nullptr */);

Finally, we measure the expectation value that we have defined earlier using the function :ref:`get_expecation_value(...)<ShellMatrix_MATMUL_get_expectation_value_A>`:

.. code-block:: cpp

    for(uint64_t s_for = 0UL; s_for < data.size(); s_for++){
        auto obs = O_matrix.get_expectation_value(states[s_for]);
        std::cout << "[main] - E_" << s_for << " = " << data[s_for] << ", O_" << s_for << " = " << obs << std::endl;
    }

The :ref:`MPI<MPI_docs>` version can be executed via:

.. code-block:: console

    mpirun -n 2 ./main -L 10 -n 5

And the :ref:`openMP<openMP_docs>` can be executed via (if activated with ``cmake -DOMP=ON ..``):

.. code-block:: console

    export OMP_NUM_THREADS=4

    ./main -L 10 -n 5


.. _Examples_SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution:

Petsc SparseMatrix **imaginary** time evolution
-----------------------------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution.rst

This :ref:`example<Examples_SparseMatrix_petsc_Hamiltonian_imaginary_time_evolution_code>` is only available if :ref:`Petsc and Slepc<Petsc_setup>` have been configured with **real** scalars (``petsc_real``) and the correct paths are set.
We are performing imaginary time evolution in a transverse field Ising model:

.. math::

   H = J \sum_{i=0}^{L-1} S^z_iS^z_{i+1} + h\sum_{i=0}^{L-1} S^x_i

.. code-block:: cpp

   Hamiltonian_NoU1 H(L);

    for(uint64_t i = 0; i < L; i++){
        H.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
        H.add_operator(h, {i}, {"Sx"});
    }


The goal is to compute the imaginary time evolution of an initial state :math:`\vert\Psi(0)\rangle` under this Hamiltonian by exploiting its sparsity:

.. math::

   \vert\Psi(\beta)\rangle = e^{-\beta (H-E_0)}\vert\Psi(0)\rangle


We need to shift the operator in the exponent such that its spectrum is semi-positive. Otherwise, the norm of the vector :math:`\vert\Psi(\beta)\rangle` grows exponentially.

We start from a random vector:

.. code-block:: cpp

	/* Generating a random initial state */
	Vec Psi;
	PetscCall(MatCreateVecs(H_matrix, &Psi, PETSC_NULLPTR));
	PetscInt start, end;
	PetscCall(VecGetOwnershipRange(Psi, &start, &end));
	std::vector<PetscInt> idx(end-start);
	std::vector<PetscScalar> data(end-start);
	std::random_device rd;
	std::mt19937 gen(rd());
	std::normal_distribution<double> dist(0, 1.);
	for(PetscInt i_for = start; i_for < end; i_for++){
		data[i_for-start] = static_cast<PetscScalar>(dist(gen));
		idx[i_for-start] = i_for;
	}
	PetscCall(VecSetValuesLocal(Psi, end-start, idx.data(), data.data(), INSERT_VALUES));
	PetscCall(VecAssemblyBegin(Psi));
	PetscCall(VecAssemblyEnd(Psi));
	data.clear();
	idx.clear();
	PetscReal norm;
	PetscCall(VecNorm(Psi,NORM_2,&norm));
	PetscCall(VecScale(Psi, 1./norm));

:ref:`Slepc<Petsc_setup>` provides the desired `matrix functions (mfn) <https://slepc.upv.es/documentation/current/docs/manualpages/MFN/index.html>`_ to perform the imaginary time evolution:

.. code-block:: cpp
	
    /* Setting up the matrix exponential with delta_beta */
    MFN mfn;
    FN f;
    PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
    PetscCall(MFNSetOperator(mfn,H_matrix));
    PetscCall(MFNGetFN(mfn,&f));
    PetscCall(FNSetType(f,FNEXP));
    PetscCall(FNSetScale(f,-delta_beta,1.0)); /* Setting the delta_beta steps */
    PetscCall(MFNSetTolerances(mfn,1e-07,PETSC_DEFAULT));
    PetscCall(MFNSetFromOptions(mfn));


The following loop performs the evolution:

.. code-block:: cpp

	for(uint64_t i = 0; i < steps; i++){

		beta += delta_beta;

		/* Time evolution PsiPrime = e^{-delta_beta (H-E0)} Psi */
		PetscCall(MFNSolve(mfn,Psi,PsiPrime));

		/* Setting Psi = PsiPrime */
		PetscCall(VecCopy(PsiPrime, Psi)); 

		/* Calculating the norm */
		PetscCall(VecNorm(Psi,NORM_2,&norm));

		/* Rescale Psi */
		PetscCall(VecScale(Psi, 1./norm));

		/* Calculating the overlap with the ground state */
		PetscCall(VecDot(Psi, Psi0, &overlap));

		/* Perform matrix-vector multiplication: PsiPrime = (H-E0) * Psi */
		PetscCall(MatMult(H_matrix, Psi, PsiPrime));

		/* Calculating the energy */
		PetscCall(VecDot(Psi, PsiPrime, &energy));

		/* Output */
		if(myrank == 0){
			std::cout << "[myrank: " << myrank << "] - [main] - beta = " << beta << ": |<Psi|Psi0>| = " << std::abs(overlap) << " - <Psi|H-E0|Psi> = " << energy << std::endl;
		}

	}

The program has to be executed with ``mpirun`` which is located in ``$MPI_DIR/bin``.
**You** can add it to ``PATH`` with ``export PATH=$MPI_DIR/bin:$PATH``.
The number of ranks are controlled with the flag ``-n 4``:

.. code-block:: console

    mpirun -n 4 ./main -L 10 -steps 10 -delta_beta 0.5 -h 0.5


.. tip::

    Alternatively, you can use the `matrix-free shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ with :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>`.

.. _Examples_SparseMatrix_petsc_Hamiltonian_real_time_evolution:

Petsc SparseMatrix **real** time evolution
------------------------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/SparseMatrix_petsc_Hamiltonian_real_time_evolution.rst

This :ref:`example<Examples_SparseMatrix_petsc_Hamiltonian_real_time_evolution_code>` is only available if :ref:`Petsc and Slepc<Petsc_setup>` have been configured with **complex** scalars (``petsc_complex``) and the correct paths are set.
We are performing real time evolution in a transverse field Ising model:

.. math::

   H = J \sum_{i=0}^{L-1} S^+_iS^-_{i+1} + S^-_iS^+_{i+1} + J_{zz} S^z_iS^z_{i+1}

The initial state is the **imbalance** state :math:`\vert Q-1;0;Q-1;0;...\rangle` which we initial as a `Petsc Vec <https://petsc.org/main/manualpages/Vec/>`_ :

.. code-block:: cpp

    /* Creating the imbalance state |Q-1;0;Q-1;0;Q-1;0;Q-1;...> as Petsc Vec */
    State state;
    for(uint64_t i_for = 0UL; i_for < L; i_for += 2){
        state.set_state(1,i_for);
    }

    /* Retrieve particle number */
    uint64_t n = state.get_particle_number(0,L);

    /* Code in between ... */

    /* Retrieving the index of state and creating the initial vector */
    PetscInt imb_index = static_cast<PetscInt>(H.get_index(state));
    PetscScalar one = 1.;
    Vec Psi;
    PetscCall(MatCreateVecs(H_matrix, &Psi, PETSC_NULLPTR));
    PetscCall(VecZeroEntries(Psi));
    PetscCall(VecSetValues(Psi, 1, &imb_index, &one, INSERT_VALUES));
    PetscCall(VecAssemblyBegin(Psi));
    PetscCall(VecAssemblyEnd(Psi));


Next, we are defining the Hamiltonian and the imbalance operator to measure its expectation value which is maximal with the just state defined above.

.. code-block:: cpp

	/* Hamiltonian */
	Hamiltonian_U1 H(L,n); 

	/* Defining XXZ Heisenberg Hamiltonian with PBC */
	for(uint64_t i = 0; i < L; i++){
		H.add_operator(1., {i,(i+1)%L}, {"S+","S-"});
		H.add_operator(1., {i,(i+1)%L}, {"S-","S+"});
		H.add_operator(Jzz, {i,(i+1)%L}, {"Sz","Sz"});
	}

	/* Observable */
	Hamiltonian_U1 O(L,n); 

	/* Imbalance operator */
	for(uint64_t i = 0; i < L; i++){
		O.add_operator(std::pow(-1.,i), {i}, {"Sz"});
	}



The goal is to compute the real time evolution of an initial state :math:`\vert\Psi(0)\rangle` under this Hamiltonian by exploiting its sparsity:

.. math::

   \vert\Psi(t)\rangle = e^{-itH}\vert\Psi(0)\rangle


Similar to the previous example, :ref:`Slepc<Petsc_setup>` provides the desired `matrix functions (mfn) <https://slepc.upv.es/documentation/current/docs/manualpages/MFN/index.html>`_ to perform the real time evolution:

.. code-block:: cpp
	
    /* Setting up the matrix exponential with delta_t */
    MFN mfn;
    FN f;
    PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
    PetscCall(MFNSetOperator(mfn,H_matrix));
    PetscCall(MFNGetFN(mfn,&f));
    PetscCall(FNSetType(f,FNEXP));
    PetscCall(FNSetScale(f,-delta_t*PETSC_i,1.0)); /* Setting the delta_t steps */
    PetscCall(MFNSetTolerances(mfn,1e-07,PETSC_DEFAULT));
    PetscCall(MFNSetFromOptions(mfn));



The following loop performs the evolution:

.. code-block:: cpp

	/* Loop of steps */
	double t = 0.;
	PetscScalar result;
	PetscReal norm;
	for(uint64_t i = 0; i < steps; i++){

		/* Perform matrix-vector multiplication: PsiPrime = O * Psi */
		MatMult(O_matrix, Psi, PsiPrime);

		/* Calculating norm and observable */
		PetscCall(VecDot(Psi, PsiPrime, &result));
		PetscCall(VecNorm(Psi,NORM_2,&norm));

		/* Time evolution PsiPrime = e^{-i delta_t H} Psi */
		PetscCall(MFNSolve(mfn,Psi,PsiPrime));

		/* Setting Psi = PsiPrime */
		PetscCall(VecCopy(PsiPrime, Psi));

		/* Output */
		if(myrank == 0){
			std::cout << "[myrank: " << myrank << "] - [main] - t = " << t << ": sqrt(<Psi|Psi>) = " << norm << " - <Psi|I|Psi> = " << result << std::endl;
		}

		t += delta_t;

	}


The program has to be executed with ``mpirun`` which is located in ``$MPI_DIR/bin``.
**You** can add it to ``PATH`` with ``export PATH=$MPI_DIR/bin:$PATH``.
The number of ranks are controlled with the flag ``-n 4``:

.. code-block:: console

    mpirun -n 4 ./main -L 10 -steps 10 -delta_t 0.5 -Jzz 2


.. tip::

    Alternatively, you can use the `matrix-free shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ with :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>`.



.. _Examples_SparseMatrix_petsc_Lindbladian:

Petsc SparseMatrix open-system dynamics
---------------------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/SparseMatrix_petsc_Lindbladian.rst

This :ref:`example<Examples_SparseMatrix_petsc_Lindbladian_code>` is only available if :ref:`Petsc and Slepc<Petsc_setup>` have been configured with **complex** scalars (``petsc_complex``) and the correct paths are set.
We are performing real time evolution in an open system described by the Lindblad master equation.
The underlying Hamiltonian is a :math:`SU(2)` symmetric Heisenberg chain with periodic boundary conditions:

.. math::

   H = \sum_{i=0}^{L-1} S^x_iS^x_{i+1} + S^y_iS^y_{i+1} + S^z_iS^z_{i+1}

Dephasing-operators :math:`S^z_i` act on all sites such that the Lindbladian is described by:

.. math::

   L[\rho] = - i [H,\rho] + \gamma \sum_{i=0}^{L-1} S^z_i\rho {S^z_i}^\dagger -\frac{1}{2}\left\{{S^z_i}^\dagger S^z_i,\rho\right\}


The model conserved the total particle number:

.. code-block:: cpp

	/* Lindbladian */
	Lindbladian_U1 Lind(L,n); 

	/* Defining the SU(2) symmetry Heisenberg model with PBC */
	double J = 1.;
	for(uint64_t i = 0; i < L; i++){
		Lind.add_operator(J, {i,(i+1)%L}, {"Sx","Sx"});
		Lind.add_operator(J, {i,(i+1)%L}, {"Sy","Sy"});
		Lind.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
		Lind.add_jump_operator(gamma, {i}, {"Sz"});
	}

	/* Hamiltonian class for the observable */
	Hamiltonian_U1 obs(L,n); 

	/* Defining the observable for measurements */
	for(uint64_t i = 0; i < L; i++){
		obs.add_operator(std::pow(-1.,i), {L-i-1}, {"Sz"});
	}



We choose the initial density matrix to be a pure state:

.. math:: 

    \rho(t=0) = \vert Q-1;0;Q-1;\dots;Q-1\rangle\langle Q-1;0;Q-1;\dots;Q-1\vert

First, we have to retrieve the correct index of :math:`\rho(t=0)` from the :ref:`Lindbladian<Operator>` using the function :ref:`get_index(state)<Operator_Functions_get_index>` that refers to the **vectorized** density matrix.
Second, we need to initial the density matrix as a `Petsc Vec <https://petsc.org/main/manualpages/Vec/>`_:

.. code-block:: cpp

	/* Creating fully imbalance density matrix rho_init = |Q-1;0;Q-1;0;...><Q-1;0;Q-1;0;...> as Petsc Vec */
	State rho_init;
	rho_init.set_range(State::template_struct::template_Q-1,0UL,2UL*L);
	for(uint64_t i_for = 0UL; i_for < L; i_for++){
		if(i_for % 2 == 1){
			rho_init.set_state(0,i_for);
			rho_init.set_state(0,i_for+L);
		}
	}

	/* Retrieve particle number */
	uint64_t n = rho_init.get_particle_number(0,L);

	/* Retrieving the index of rho_init to create the initial vector */
	PetscInt imb_index = static_cast<PetscInt>(Lind.get_index(rho_init));
	PetscScalar one = 1.;
	Vec rho;
	PetscCall(MatCreateVecs(Lind_matrix, &rho, PETSC_NULLPTR));
	PetscCall(VecZeroEntries(rho));
	PetscCall(VecSetValues(rho, 1, &imb_index, &one, INSERT_VALUES));
	PetscCall(VecAssemblyBegin(rho));
	PetscCall(VecAssemblyEnd(rho));

Besides defining the Lindbladian, we define the imbalance operator is given by:

.. math::

   O = \sum_i (-1)^{i} S_i^z

.. code-block:: cpp

	/* Lindbladian */
	Lindbladian_U1 Lind(L,n); 

	/* Defining the SU(2) symmetry Heisenberg model with PBC */
	double J = 1.;
	for(uint64_t i = 0; i < L; i++){
		Lind.add_operator(J, {i,(i+1)%L}, {"Sx","Sx"});
		Lind.add_operator(J, {i,(i+1)%L}, {"Sy","Sy"});
		Lind.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
		Lind.add_jump_operator(gamma, {i}, {"Sz"});
	}

	/* Hamiltonian class for the observable */
	Hamiltonian_U1 obs(L,n); 

	/* Defining the observable for measurements */
	for(uint64_t i = 0; i < L; i++){
		obs.add_operator(std::pow(-1.,i), {L-i-1}, {"Sz"});
	}


We define a `matrix functions (mfn) <https://slepc.upv.es/documentation/current/docs/manualpages/MFN/index.html>`_, similar to the previous example, to perform the time evolution via the matrix exponential applied to the initial density matrix :math:`e^{tL}\vert \rho(t=0)\rangle`:


.. code-block:: cpp

	/* Setting up the matrix exponential with deltaT */
	MFN mfn;
	FN f;
	PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
	PetscCall(MFNSetOperator(mfn,Lind_matrix));
	PetscCall(MFNGetFN(mfn,&f));
	PetscCall(FNSetType(f,FNEXP));
	PetscCall(FNSetScale(f,deltaT,1.0)); /* Setting the deltaT steps */
	PetscCall(MFNSetTolerances(mfn,1e-8,PETSC_DEFAULT));
	PetscCall(MFNSetFromOptions(mfn));

Lastly, we put everything (time-evolution, entropy, and imbalance measurement) together in a simple for loop.
The imbalance is computed using the function :ref:`trace_over_density_matrix(rho)<Operator_Functions_Miscellaneous_trace_over_density_matrix>` and the von Neumann entropy is computed with :ref:`von_Neumann_entropy_from_PetscVec(rho)<Functions_von_Neumann_entropy_from_PetscVec>`:

.. code-block:: cpp


	/* Loop of steps */
	double t = 0.;
	for(uint64_t i = 0; i < steps; i++){
		t += delta_t;

		/* Time evolution rhoPrime = e^{delta_t L} rho */
		PetscCall(MFNSolve(mfn,rho,rhoPrime));

		/* Setting rho = rhoPrime */
		PetscCall(VecCopy(rhoPrime, rho));

		/* Retrieving the density matrix as a Eigen matrix on rank 0 */
		double S = von_Neumann_entropy_from_PetscVec(rho,0);

		/* Retrieving the imbalance from the operator obs */
		PetscScalar imbalance = obs.trace_over_density_matrix(rho);

		/* Output */
		if(myrank == 0){
			std::cout << "[myrank: " << myrank << "] - [main] - t = " << t << ": S/Smax = " << S/std::log(dim)*2. << ", - imbalance = " << imbalance << std::endl;
		}
	}

The program has to be executed with ``mpirun`` which is located in ``$MPI_DIR/bin``.
**You** can add it to ``PATH`` with ``export PATH=$MPI_DIR/bin:$PATH``.
The number of ranks are controlled with the flag ``-n 4``:

.. code-block:: console

    mpirun -n 4 ./main -L 10 -steps 20 -delta_t 0.5 -gamma 0.5


.. tip::

    Alternatively, you can use the `matrix-free shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ with :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>`.

.. _Examples_ShellMatrix_petsc_Hamiltonian:

Petsc ShellMatrix ground state
------------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/ShellMatrix_petsc_Hamiltonian.rst


We determine the ground state energy using the :ref:`matrix-free shell<Operator_Functions_Frontends_create_PetscShellMatrix>`.
The :ref:`example<Examples_ShellMatrix_petsc_Hamiltonian_code>` is only available if :ref:`Petsc and Slepc<Petsc_setup>` have been configured and ``petsc_real`` is available.

We focus on an one-dimensional XXZ-Heisenberg model with periodic boundary conditions:

.. math::

   H = \sum_{i=0}^{L-1} \frac{1}{2}\left(S^+_iS^-_{i+1} +S^-_iS^+_{i+1}\right) + J_z S^z_iS^z_{i+1}


.. code-block:: cpp

    Hamiltonian_U1 H(L,n);

    for(uint64_t i = 0; i < L; i++){
        H.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
        H.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
        H.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
    }

Next, we retrieve the `matrix shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ with :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>`:


.. code-block:: cpp

    auto H_matrix = H.create_PetscShellMatrix();


The `matrix shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ can be used for all :ref:`Slepc routines<Petsc_setup>` that are based on matrix multiplications such as `Krylov space methods <https://en.wikipedia.org/wiki/Krylov_subspace>`_.
This includes the `Eigenvalue Problem Solver (eps) <https://slepc.upv.es/documentation/current/docs/manualpages/EPS/index.html>`_ to determine the ground state:

.. code-block:: cpp

    /* Setting up the eigenvalue problem */
    EPS eps;
    PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
    PetscCall(EPSSetOperators(eps, H_matrix, NULL));
    PetscCall(EPSSetProblemType(eps, EPS_HEP));
    PetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
    PetscCall(EPSSetDimensions(eps, 1, PETSC_DEFAULT, PETSC_DEFAULT));

    /* Calculate ground state energy */
    PetscCall(EPSSolve(eps));
    PetscInt nconv;
    PetscCall(EPSGetConverged(eps, &nconv));
    PetscScalar E0, E0_imag;
    PetscCall(EPSGetEigenvalue(eps,0, &E0, &E0_imag));



The program has to be executed with ``mpirun`` which is located in ``$MPI_DIR/bin``.
**You** can add it ``PATH`` with ``export PATH=$MPI_DIR/bin:$PATH``.
The number of ranks are controlled with the flag ``-n 4``:

.. code-block:: console

    mpirun -n 4 ./main -L 10 -n 5 -Jz 0.5



.. tip::

    Alternatively, you store the sparse matrix in a `MATMPIAIJ <https://petsc.org/main/manualpages/Mat/MATMPIAIJ/>`_ format with :ref:`create_PetscSparseMatrix()<Operator_Functions_Frontends_create_PetscSparseMatrix>`.


.. _Examples_ShellMatrix_petsc_Lindbladian:

Petsc ShellMatrix steady state
------------------------------

.. toctree::
    :maxdepth: 1
    :titlesonly:
    
    Examples/ShellMatrix_petsc_Lindbladian.rst


We determine the steady state and decay rate using the :ref:`matrix-free shell<Operator_Functions_Frontends_create_PetscShellMatrix>` for a Lindbladian with non-hermitian jump operators.
The :ref:`example<Examples_ShellMatrix_petsc_Lindbladian_code>` requires the complex version (``petsc_complex``) of :ref:`Petsc and Slepc<Petsc_setup>`.

The Hamiltonian is described by an one-dimensional XXZ-Heisenberg model with periodic boundary conditions:

.. math::

   H = \sum_{i=0}^{L-1} \frac{1}{2}\left(S^+_iS^-_{i+1} +S^-_iS^+_{i+1}\right) + J_z S^z_iS^z_{i+1}

and we define two jump operators acting on the first site:

.. math::

   L[\rho] = - i [H,\rho] + \gamma \sum_{L\in\{S^+_0 ,S^-_0\}} L\rho L^\dagger -\frac{1}{2}\left\{L^\dagger L,\rho\right\}


Note that the Lindbladian does not conserve the particle number and dimension of the Lindbladian is :math:`Q^{2L}`.

.. code-block:: cpp

	/* Operator class */
	Lindbladian_NoU1 Lind(L); 

	/* Simply XXZ-Heisenberg model with periodic boundary conditions */
	for(uint64_t i = 0; i < L; i++){
		Lind.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
		Lind.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
		Lind.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
	}

	/* Non-hermitian jump operators */
	Lind.add_jump_operator(gamma, {0}, {"S+"});
	Lind.add_jump_operator(gamma, {0}, {"S-"});



Similar to the previous example, we retrieve the `matrix shell <https://petsc.org/release/manualpages/Mat/MATSHELL/>`_ with :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>` and set up the eigenvalue problem:


.. code-block:: cpp

	/* Creating the matrix shell for matrix-free for multiplication */
	auto Lind_matrix = Lind.create_PetscShellMatrix();
	
	/* Setting up the eigen problem solver */
	EPS eps;
	PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
	PetscCall(EPSSetOperators(eps, Lind_matrix, NULL));
	PetscCall(EPSSetProblemType(eps, EPS_NHEP));
	PetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL));
	PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
	PetscCall(EPSSetDimensions(eps, nev, ncv, PETSC_DEFAULT));

	/* Define tolerance, last argument refers to the maximal number of iterations */
	PetscCall(EPSSetTolerances(eps,tol,10000));

	/* Setting the monitoring function */
	PetscCall(EPSMonitorSet(eps,&monitor,NULL,NULL));

	/* Solve the system */
	PetscCall(EPSSolve(eps));




Eigenvalues of the Lindbladian are complex with a real part smaller or equal to 0.
Crucially, Krylov-space method converges much better towards the ones with the largest magnitude.
However, steady states and the smallest decay rates (eigenvalues with the largest real part) have the smallest magnitude, making the convergence difficult in some cases.
Hence, depending on the model, changing some parameters of the algorithm helps enhance the convergence.
One such trick employed here is a spectral transform by defining a shift on the real axis.
Also, changing the number of stored basis vectors for the Krylov space (**ncv**)  can help.
The following things can be changed at runtime:

* **L**: system size

* **gamma**: coupling strength

* **Jz**: interaction strength

* **nev**: number of eigenvalues

* **tol**: tolerance for convergence criteria

* **ncv**: stored basis vectors

* **shift**: spectral shift

Note that the program has to be executed with ``mpirun`` which is located in ``$MPI_DIR/bin``.
**You** can add it ``PATH`` with ``export PATH=$MPI_DIR/bin:$PATH``.
The number of ranks are controlled with the flag ``-n 4``:

.. code-block:: console

    mpirun -n 4 ./main -L 5 -gamma 1 -nev 3


.. tip::

    Alternatively, you store the sparse matrix in a `MATMPIAIJ <https://petsc.org/main/manualpages/Mat/MATMPIAIJ/>`_ format with :ref:`create_PetscSparseMatrix()<Operator_Functions_Frontends_create_PetscSparseMatrix>`.