Examples

We illustrate the usage using various examples that include the different Operator types (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 CMake to build the source files:

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:

State and Basis class

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

/* 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 Basis class with the Basis Iterator. The Basis class is built on top of the Container class that built on top of the State class. It is constructed by defining the system size L, the particle sector n, and the number of subsystems Npart. The Iterator class, it below, allows an easy handling of the Basis:

/* 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;

}

Operator class

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

/* 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>;

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

\[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}\]
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 create_DenseMatrix() or in a sparse format, see SparseMatrix, with create_SparseMatrix().

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

Both functions are compatible with openMP. This can be activated with cmake -DOMP=ON ...

Eigen full spectrum

Functions in this example are available if the EIGEN_DIR is set correctly. A detailed guide to set up Eigen can be found here. Eigen should be used for dense computations such as determining the full spectrum of an operator. In this example, we look at an one-dimensional XXX-Heisenberg model with periodic boundary conditions and an alternating field in \(S^z_i\).

\[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:

\[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\}\]
/* 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 Eigen matrix can be retrieved with create_EigenDenseMatrix(). The functions is compatible with openMP.

auto Lind_matrix = Lind.create_EigenDenseMatrix();

The returned matrix can be used in all Eigen routines. The following code determines all eigen values:

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

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

auto eigenvalues = solver.eigenvalues();

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

export OMP_NUM_THREADS=4
export MKL_NUM_THREADS=4

./main -L 10 -n 5

ShellMatrix ground state

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

Given an Operator, the ShellMatrix provides a matrix-vector multiplication (without storing the matrix) that can be used to generate Krylov subspaces to compute ground states and other properties. States are represents by the Vector class that – as well – supports MPI and openMP.

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

\[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\]
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 ShellMatrix is simply generated with create_ShellMatrix():

Hamiltonian_U1 O(L,n);

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

auto O_matrix = O.create_ShellMatrix();

We have implemented the Lanczos algorithm to compute the ground state and excitations with ShellMatrix class. It is simply called by passing the Hamiltonian to the function lanczos(H,…).

/* 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 get_expecation_value(…):

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 MPI version can be executed via:

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

And the openMP can be executed via (if activated with cmake -DOMP=ON ..):

export OMP_NUM_THREADS=4

./main -L 10 -n 5

Petsc SparseMatrix imaginary time evolution

This example is only available if Petsc and Slepc 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:

\[H = J \sum_{i=0}^{L-1} S^z_iS^z_{i+1} + h\sum_{i=0}^{L-1} S^x_i\]
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 \(\vert\Psi(0)\rangle\) under this Hamiltonian by exploiting its sparsity:

\[\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 \(\vert\Psi(\beta)\rangle\) grows exponentially.

We start from a random vector:

/* 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));

Slepc provides the desired matrix functions (mfn) to perform the imaginary time evolution:

/* 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:

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:

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

Tip

Alternatively, you can use the matrix-free shell with create_PetscShellMatrix().

Petsc SparseMatrix real time evolution

This example is only available if Petsc and Slepc 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:

\[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 \(\vert Q-1;0;Q-1;0;...\rangle\) which we initial as a Petsc Vec :

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

/* 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 \(\vert\Psi(0)\rangle\) under this Hamiltonian by exploiting its sparsity:

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

Similar to the previous example, Slepc provides the desired matrix functions (mfn) to perform the real time evolution:

/* 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:

/* 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:

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

Tip

Alternatively, you can use the matrix-free shell with create_PetscShellMatrix().

Petsc SparseMatrix open-system dynamics

This example is only available if Petsc and Slepc 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 \(SU(2)\) symmetric Heisenberg chain with periodic boundary conditions:

\[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 \(S^z_i\) act on all sites such that the Lindbladian is described by:

\[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:

/* 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:

\[\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 \(\rho(t=0)\) from the Lindbladian using the function get_index(state) that refers to the vectorized density matrix. Second, we need to initial the density matrix as a Petsc Vec:

/* 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:

\[O = \sum_i (-1)^{i} S_i^z\]
/* 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), similar to the previous example, to perform the time evolution via the matrix exponential applied to the initial density matrix \(e^{tL}\vert \rho(t=0)\rangle\):

/* 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 trace_over_density_matrix(rho) and the von Neumann entropy is computed with von_Neumann_entropy_from_PetscVec(rho):

/* 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:

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

Tip

Alternatively, you can use the matrix-free shell with create_PetscShellMatrix().

Petsc ShellMatrix ground state

We determine the ground state energy using the matrix-free shell. The example is only available if Petsc and Slepc have been configured and petsc_real is available.

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

\[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}\]
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 with create_PetscShellMatrix():

auto H_matrix = H.create_PetscShellMatrix();

The matrix shell can be used for all Slepc routines that are based on matrix multiplications such as Krylov space methods. This includes the Eigenvalue Problem Solver (eps) to determine the ground state:

/* 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:

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

Tip

Alternatively, you store the sparse matrix in a MATMPIAIJ format with create_PetscSparseMatrix().

Petsc ShellMatrix steady state

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

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

\[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:

\[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 \(Q^{2L}\).

/* 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 with create_PetscShellMatrix() and set up the eigenvalue problem:

/* 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:

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

Tip

Alternatively, you store the sparse matrix in a MATMPIAIJ format with create_PetscSparseMatrix().