Operator

Summary

The class can be built from the State class and the Basis class with number conservation.
The interface is similar to the ITensor package and allows a straightforward usage.
It provides different formats to extract the matrix such as Eigen and Petsc.
It handles Hamiltonian and Lindbladian systems

Operator Type

The operator can be defined on space of wavefunction and on the space of the density matrices. The latter is required when we deal with dissipative system and the dynamics is described the Lindblad formalism.

Both can be used with and without number conservation and all four cases are predefined in the header file DanceQ.h (L, n refer to the system size, particle number – if desired):

danceq::Lindbladian_U1(L,n)

danceq::Lindbladian_NoU1(L)

danceq::Hamiltonian_U1(L,n)

danceq::Hamiltonian_NoU1(L)

The class Operator class has three template parameters (Basis, ScalarType, and OpertorType). The OperatorType is a enum that distinguishes the two different cases:

enum danceq::OperatorType

Distinguishes between Hamiltonian and Lindbladian systems.

The Lindbladian type acts on the space of vectorized density matrices:

\[\rho_{\text{vectorized}}[i\cdot d +j] = \rho_{\text{matrix}}[i,j]\]

Here, \(d\) is the underlying Hilbert space dimension, \([X]\) refers to the index of the vectorized density matrix, \([X,Y]\) refer to the row and column entries of the density matrix in matrix form. A Hamiltonian operator can be used to compute the trace over a vectorized density matrix with the function trace_over_density_matrix(rho_vectorized).

Values:

enumerator Hamiltonian: Defines a Hamiltonian operator acting on the space of wavefunctions.

enumerator Lindbladian: Defines a Lindbladian operator acting on the space of density matrices.

The Lindbladian acts on the space of density matrices \(\rho\):

\[\frac{\text{d}\rho}{\text{d}t} = -i[H,\rho] + \sum_{l} l\rho l^\dagger - \gamma\frac{1}{2}l^\dagger l\rho - \gamma\frac{1}{2}\rho l^\dagger l\]

Basis

The underlying basis of the operator can be set to the State class with dimension \(\mathrm{Q}^L\) or to the Basis class with dimension from Eq. (2):

/* Operator using the Basis class, system size L, particle number n */
using MyOperatorU1=danceq::internal::Operator<Basis,double,Hamiltonian> Operator;
MyOperatorU1 H_U1(L,n);

/* Operator using the State class, system size L */
using MyOperatorNoU1=danceq::internal::Operator<State,double,Hamiltonian>;
MyOperatorNoU1 H_NoU1(L);

The second template parameter defines the scalar, e.g., complex and real.

Warning

MaxSites defined in the State class has to be larger or equal the system size if we deal with Hamiltonian systems and at least double the size if we deal Lindbladian systems.

Input

The interface is similar to the ITensor and provides an easy way to create arbitrary linear operators from tensor products with the add_operator(…) function. The following example creates a one-dimensional Heisenberg model of size L with n particles, coupling constant J, a local alternating field h, and a dissipation rate gamma:

/* system size */
uint64_t L = 10;

/* particle number */
uint64_t n = 5;

/* field strength */
double h = 0.1;

/* coupling strength */
double J = 1.;

/* dissipation rate */
double gamma = 0.01;

/* defining the Lindbladian with number conservation */
Lindbladian_U1(L,n);

/* adding the Hamiltonian terms */
for(uint64_t i = 0; i < L; i++){
    L.add_operator(J/2., {i,(i+1)%L}, {"S+","S-"});
    L.add_operator(J/2., {i,(i+1)%L}, {"S-","S+"});
    L.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
    L.add_operator(std::pow(-h,i), {i}, {"Sz"});
}

/* jump operator acting on site 0 */
L.add_jump_operator(gamma, {0}, {"Sz"});

A list of default operators is already defined. New local operators of size \(\textbf{Q}\times\textbf{Q}\) acting on a single site can be defined using define_operator(…) with a std::array< std::array<ScalarType, Q>, Q>.

Matrix formats

The class comes with different formats to represent the operator as a matrix using the provided Basis. Besides native cpp containers, frontends for Eigen and Petsc are available if they are included by CMake

The dense matrix can be retrieved using a two-dimensional std::vector<std::vector<ScalarType>>. We further provided an advanced SparseMatrix calss that can be used for simple routines such as Matrix-Vector multiplication. The class works together with the provided Vector class that is a primitive extension of std::vector. It is capable of primitive routines and provides distributed memory if MPI is available. Together with the SparseMatrix calss naive Krylov space methods can be used such as ground state calculation. No additional third-party libraries like Eigen or Petsc and Slepc are required. For more complex operators, we recommend working with more advanced linear algebra packages as listed below:

Eigen can be used for dense applications (e.g., calculating the full spectrum) and provides access to LAPACK libraries such as Intel’s MKL. For sparse operations (e.g., Krylov space methods), Petsc and Slepc provide a powerful MPI-based approach allowing scalable computation across multiple nodes! Additionally, we have added a matrix-free implementation with Petsc to reduce the memory capacity. The application does not store the matrix and computes its elements on-the-fly.

Available matrix formats

create_DenseMatrix(): Dense matrix using two-dimensional std::vector
create_EigenDenseMatrix(): Dense matrix using Eigen
create_SparseMatrix(): Sparse matrix using the SparseMatrix class
create_ShellMatrix(): Matrix-free multiplication using the ShellMatrix class
create_PetscSparseMatrix(): Sparse matrix using Petsc
create_PetscShellMatrix(): Matrix-free multiplication using Petsc with the ShellMatrix class

Application to a product states

The core function which is used to create any matrix format is apply(…). It returns the elements (indices and scalars) of i th column when applied to in the i th product state of the underlying basis.

In order to optimize the computation, we exploit a similar divide-and-conquer ansatz acting on non-local subclusters. When a new operator string is added by add_operator(…) it is only stored in operator_list and is not directly used for the computation. The apply(…) function (and therefore all other more complex functions) is working with the data stored in sparse_matrices_for_computation which is created by the merge_terms(…). It bundles all operator strings from operator_list in sparse_matrices_for_computation. You do not need to call merge_terms(…) since it is called every time a new operator string was added by add_operator(…).

The merge_terms(…) function identifies clusters of maximal size Ntensor_max and integrates all operator strings that are fully supported on such a cluster. Their additive contribution is stored in a sparse format – similar to the SparseMatrix – within sparse_matrices_for_computation. The default value of Ntensor_max is chosen such that the dimension of the underlying sparse matrix does not exceed 2048. Ntensor_max can be set by calling set_Ntensor_max(…).

Full example

 #include "Operator.h"

 /* Defining the state class with a maximal length of 256 sites and a local dimension of 4 */
 using State=danceq::internal::State<256,4>;

 /* Defining the Basis class from the State class */
 using Basis=danceq::internal::BasisU1<State>;

 /* Defining the Operator from the Basis class */
 using MyOperator=danceq::internal::Operator<Basis,double,Hamiltonian>;

 /* Defining the sparse matrix class using the scalar from the Operator class */
 using SparseMatrix=danceq::internal::SparseMatrix<Operator::scalartype>;


 int main(int argc, char *argv[]) {

     /* Input parameter */
     uint64_t L = 10;
     uint64_t n = 5;
     double h = 0.1;
     double J = 1.;

     /* Operator class */
     MyOperator H(L,n);

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

     /* Simply XXX-Heisenberg model with a local field along the z-axis */
     for(uint64_t i = 0; i < L; i++){
         H.add_operator(J/2, {i,(i+1)%L}, {"S+","S-"});
         H.add_operator(J/2, {i,(i+1)%L}, {"S-","S+"});
         H.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
         H.add_operator(std::pow(-h,i), {i}, {"Sz"});
     }

     /* Dense matrix using std::vector<std::vector<Operator::scalartype>> */
     auto H_dense = H.create_DenseMatrix();

     /* Sparse matrix using Sprase_Matrix class */
     auto H_sparse = H.create_SparseMatrix();

     /* Dense matrix with Eigen using the type Operator::EigenMatrixType */
     auto H_eigen = H.create_EigenDenseMatrix();

     /* Sparse matrix with Petsc Mat type MATMPIAIJ */
     auto H_petsc = H.create_PetscSparseMatrix();

     /* Sparse matrix with Petsc using a matrix-free shell */
     auto H_petsc_shell = H.create_PetscShellMatrix();

     return 0;
 }

Content

Detailed documentation for the class can be found in Data Members, Functions, and Autodocs Operator.